Glacial AMOC shoaling despite vigorous tidal dissipation: vertical stratification matters

During the Last Glacial Maximum (LGM), tidal dissipation was about 3-fold higher than today, which could have led to a considerable increase in vertical mixing. This increase might have enhanced the glacial Atlantic Meridional Overturning Circulation (AMOC), contradicting the shoaled AMOC indicated by paleoproxies. Here, we conduct ocean model simulations to investigate the impact of background climate conditions and tidal mixing on the AMOC during the LGM. We successfully reproduce the stratified ocean characteristics of the LGM by accurately simulating the elevated salinity of the deep sea and the rapid temperature decrease in the ocean's upper layers. Our findings indicate that the shoaled glacial AMOC is mainly due to strong glacial-ocean stratification, regardless of enhanced tidal dissipation. However, glacial tidal dissipation plays a critical role in the intensification of Antarctic Bottom Water (AABW) during the LGM. Given the critical role of the AMOC in (de-)glacial climate evolution, our results highlight the complex interactions of ocean stratification and tidal dissipation that have been neglected so far.

1 Introduction

The Atlantic Meridional Overturning Circulation (AMOC) transports heat over large distances and is therefore an essential component of the Earth's present and past climate systems (Ganopolski and Rahmstorf, 2001; Gordon, 1986; Rahmstorf, 1996; Stute et al., 2001). A major focus of paleoceanography involves understanding the contribution of the AMOC to glacial–interglacial climate change (Boyle and Keigwin, 1987; Broecker and Hemming, 2001; Clark et al., 2002; Knorr and Lohmann, 2003; Knorr et al., 2021).

The state of deep-water formation during the Last Glacial Maximum (LGM) has been discussed in the paleoclimate literature. Based on water mass properties, North Atlantic Deep Water (NADW) formation was shallower (Butzin et al., 2005; Curry and Oppo, 2005; Duplessy et al., 1988; Ferrari et al., 2014; Hesse et al., 2011; Lippold et al., 2012; Lund et al., 2011; Lynch-Stieglitz et al., 2007; Muglia et al., 2018; Skinner et al., 2017) but not much weaker (McManus et al., 2004; Sarnthein et al., 1994). At the same time, Antarctic Bottom Water (AABW) export from the Southern Ocean increased (Ledbetter and Johnson, 1976; Negre et al., 2010; Robinson et al., 2005). The salinity of glacial AABW may have been much greater than that observed today, leading to enhanced stratification of the glacial ocean between the upper and lower cells (Adkins et al., 2002; Bouttes et al., 2009; Francois et al., 1997; Jansen, 2017; Klockmann et al., 2016; Knorr et al., 2021; Lund et al., 2011; Stein et al., 2020; Watson and Garabato, 2006). Several modeling studies suggest a physical basis for the shoaled glacial AMOC, likely caused by changes in Southern Ocean sea ice (Baker et al., 2020; Butzin et al., 2005; Ferrari et al., 2014; Jansen and Nadeau, 2016; Marzocchi and Jansen, 2017; Nadeau et al., 2019; Sun et al., 2018, 2020; Watson et al., 2015) or terrestrial-ice input (Miller et al., 2012).

Coupled ocean–atmosphere model simulations of the LGM climate reveal a broad spectrum of results, showing considerable disagreement regarding whether the AMOC was weaker or stronger compared to present-day (PD) conditions (Kageyama et al., 2021; Knorr et al., 2021; Otto-Bliesner et al., 2007; Weber et al., 2007; Zhang et al., 2013). However, a critical factor often overlooked in these analyses is the significantly enhanced tidal dissipation that occurred during the LGM (Arbic et al., 2004b; Egbert et al., 2004; Green, 2010; Griffiths and Peltier, 2008, 2009; Wilmes and Green, 2014). Incorporating this element into models, as demonstrated in the research by Schmittner et al. (2015) and Wilmes et al. (2019), leads to a notable finding: both the depth and the strength of the AMOC during the LGM are substantially increased when changes in tidal dissipation are taken into account. This suggests a pivotal role for tidal mixing in shaping the LGM's ocean circulation dynamics.

Currently, tides provide about half, or 1 TW (1 TW = 10¹² W), of the energy required to maintain the global meridional overturning circulation (MOC) (Ferrari and Wunsch, 2009; Wunsch and Ferrari, 2004). Numerous studies have suggested a significant intensification of tides due to a 120–130 m drop in global mean sea levels and the exposure of continental shelves during the LGM (Arbic et al., 2004b; Egbert et al., 2004; Green, 2010; Griffiths and Peltier, 2008, 2009; Wilmes and Green, 2014). This exposure reduces effective damping, leading to an increase in tides. Additionally, there was greater tidal dissipation in the deep-ocean interior than on the continental shelves during the LGM. This amplified tidal dissipation may have been a critical factor in driving a more vigorous glacial AMOC compared to current levels, as postulated by Green et al. (2009), Schmittner et al. (2015), and Wilmes et al. (2019). Therefore, changes in tidal dissipation play an important role and should not be neglected in paleoclimate simulations (Schmittner et al., 2015). It is noteworthy that Wilmes et al. (2021) achieved a relatively shoaled AMOC through the artificial reduction of meridional moisture flux and precipitation at high latitudes. However, to date, no research has directly demonstrated a shoaled AMOC under realistic LGM forcing conditions, despite the presence of enhanced glacial tidal dissipation.

This study has three primary objectives:

1. to reproduce a stratified glacial ocean and a shoaled AMOC under actual LGM forcing conditions, considering both increased local glacial dissipation and far-field tidal dissipation;

2. to analyze the reasons for a shoaled AMOC despite the presence of enhanced tidal dissipation during the LGM;

3. to compare modeled ocean circulation with paleoclimate reconstructions.

To achieve these goals, we employ a global OGCM (ocean general circulation model) to generate a series of ocean circulation scenarios. These scenarios are driven by both LGM and PD surface forcing, as well as varying degrees of tidal mixing. Previous studies have already explored the role of increased glacial-ocean stratification in causing a shallower AMOC during the LGM (Jansen and Nadeau, 2016; Jansen, 2017). In this context, our analysis underscores that, despite the nearly 3-fold intensification of tidal dissipation during the LGM, enhanced stratification still plays a dominant role in maintaining a shoaled glacial AMOC.

2 Materials and methods

2.1 Tidal model

The global tidal model is based on the Finite Volume Community Ocean Model (FVCOM), which uses an unstructured finite-volume model with triangular meshes (Chen et al., 2003). The tidal model solves the following equation:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathbf{U}}{\partial t}}& +f\times \mathbf{U}+\mathbf{U}\cdot \mathrm{\nabla }\mathbf{u}=-gH\mathrm{\nabla }\left(\mathit{\zeta }-\mathit{\alpha }{\mathit{\zeta }}_{\mathrm{EQ}}-{\mathit{\zeta }}_{\mathrm{SAL}}\right)\\ & -a_H{\mathrm{\nabla }}^{\mathrm{2}}\mathbf{U}+{\mathbf{D}}_{\mathrm{BL}}+{\mathbf{D}}_{\mathrm{IT}},\end{array}\end{array}$$

where u is the horizontal velocity; U=uH is the horizontal transport speed; f is the Coriolis parameter; g is the gravitational acceleration; ζ is the instantaneous tide level; ζ_EQ is the equilibrium tide level (Hendershott, 1972); α is the body tide Love number; ζ_SAL is the term for gravitational self-attraction and loading tides, implemented using an iterative method (Arbic et al., 2004a; Egbert et al., 2004); and A_H is the horizontal turbulent-eddy-viscosity coefficient. Momentum is dissipated through two processes. First, we use a bottom friction term that is quadratic in velocity:

$$\begin{array}{}\text{(2)}& {\mathbf{D}}_{\mathrm{BL}}=-C_d\mathbf{u}\left|\mathbf{u}\right|,\end{array}$$

where the bottom friction coefficient (C_d) is set at 0.0025. Second, we use D_IT as the internal wave drag, i.e., the linear transfer of energy to internal waves, based on Zaron and Egbert (2006). It is expressed as follows:

$$\begin{array}{}\text{(3)}& {\mathbf{D}}_{\mathrm{IT}}=\mathrm{\Gamma }H(\mathrm{\nabla }H{)}^{\mathrm{2}}{\displaystyle \frac{N_{\mathrm{b}}\overline{N}}{\mathrm{8}{\mathit{\pi }}^{\mathrm{2}}\mathit{\omega }}}\mathbf{u},\end{array}$$

where Γ= 50 is the scaling factor. N_b and $\overline{N}$ are the buoyancy frequency at the seafloor and the depth-averaged vertical value, respectively, and are both derived from our simulations using version 2 of the Finite-volumE Sea ice–Ocean Model (FESOM2.0; see description in Sect. 2.2). It is noteworthy that the tidal dissipation obtained from the tidal model may further influence the N² value obtained from FESOM2.0, leading to certain sensitivities. We employ an iterative process to eliminate these sensitivities, with the detailed iterative process provided in Appendix A1. Moreover, ω is the tidal frequency of the M2 tide, a major tidal constituent. The tidal model utilizes four major tidal constituents (M2, S2, K1, and O1), accounting for more than 94 % of today's dissipation (Egbert and Ray, 2003). The experiments are executed for a total of 30 d, with the final 20 d used for harmonic analysis.

The resolution of the model ranges from 10 to 40 km, with a higher resolution for shallow waters and areas with significant water depth changes. For the PD triangular mesh, there are 422 932 nodes and 817 641 cells. For the LGM, the numbers are 323 101 and 624 641, respectively. The term “node” refers to each vertex of the unstructured triangular mesh, whereas “cell” denotes each triangle formed by connecting these nodes.

Here, we calculate the bottom friction dissipation (D_BL) and the internal-tide dissipation (D_IT) due to the linear transfer of energy to internal waves:

$$\begin{array}{}\text{(4)}& D_{\mathrm{IT}}=<{\mathit{\rho }}_{\mathrm{0}}\mathbf{u}\cdot {\mathbf{D}}_{\mathrm{IT}}>,D_{\mathrm{BL}}=<{\mathit{\rho }}_{\mathrm{0}}\mathbf{u}\cdot {\mathbf{D}}_{\mathrm{BL}}>,\end{array}$$

where the reference density is set to 1035 kg m⁻³ and the angle brackets (“<” and “>”) denote the tide period of the respective tidal constituent.

Figure 1 presents the bottom friction dissipation (D_BL) and internal-tide dissipation (D_IT) during the PD and the LGM. A significant difference in their distributions can be observed: bottom friction dissipation (D_BL) is mainly concentrated in shallow sea areas, while internal-tide dissipation (D_IT) is primarily found in deep-sea regions with notable topographic variations, such as mid-ocean ridges. Moreover, compared to the PD, primary regions of energy dissipation shifted from shallow to deep seas during the LGM. In the PD, dissipation below 500 m depth accounts for 1.86 TW, which is 62 % of the total dissipation. In contrast, during the LGM, dissipation below 500 m amounts to 1.35 TW, accounting for only 29 % of the total dissipation (Table 1). These data are consistent with previous research findings (Arbic et al., 2004a; Egbert et al., 2004; Green, 2010; Griffiths and Peltier, 2008, 2009; Wilmes and Green, 2014).

Figure 1 Global distributions of bottom friction dissipation (D_BL) and internal-tide dissipation (D_IT) during the PD and LGM.

Table 1 Global and sub-500 m distributions of D_IT and D_BL during the PD and LGM (measured in terawatts).

2.1.1 Bathymetry

The PD bathymetry data comes from the 1 min RTopo-2 database (Schaffer et al., 2016). For the LGM bathymetry, we use sea-level data from version 1.2 of ICE-5G (VM2 L90) (Peltier, 2004). Notably, tidal dissipation derived using ICE-6G is actually weaker than that obtained using ICE-5G (Wilmes et al., 2019; Wilmes et al., 2021). Here, we have chosen to use ICE-5G to investigate whether the AMOC during the LGM would have been affected by these stronger tidal conditions. The sea-level difference between the PD and LGM is calculated by subtracting the PD sea levels from the LGM sea levels obtained from the respective ICE-5G dataset. The low-resolution paleo-sea-level changes (1° horizontal resolution) are then interpolated to the RTopo-2 grid and added to the PD RTopo-2 bathymetry in order to retain the high-resolution topographic features. Finally, we interpolate the high-resolution bathymetry to the unstructured triangular mesh of the tidal model.

2.1.2 Tidal-model validation

The harmonically analyzed amplitudes (complex sinusoids) are used to evaluate the elevations. Simulated sinusoids ($\widehat{\mathit{\zeta }}$) were interpolated to the 0.17° grid of TPXO9.v1 and compared to the reference sinusoids ($\widehat{{\mathit{\zeta }}_R}$) by evaluating the spatially averaged root-mean-square (RMS) error (Δζ),

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }\mathit{\zeta }=\sqrt{{\displaystyle \frac{\iint \mathrm{|}\widehat{\mathit{\zeta }}-\widehat{{\mathit{\zeta }}_R}{\mathrm{|}}^{\mathrm{2}}\mathrm{d}A}{\mathrm{2}\iint \mathrm{d}A}}}.\end{array}$$

The RMS errors were calculated for four tidal constituents and determined separately for deep-water regions (depths > 500 m) and shallow shelf seas. The results are presented in Table 2. Meanwhile, Table 3 compares the M2 RMS error with the RMS errors of other forward tidal models.

Table 2 Tidal-model RMS errors for four tidal constituents.

Table 3 Comparison of M2 RMS errors in forward tidal models. Note that NA stands for not available.

^* Note that differences in the models selected across different studies may influence the results.

2.2 Ocean model

FESOM2.0 (Danilov et al., 2017), the ocean component of the Alfred Wegener Institute (AWI) Earth System Model (Sidorenko et al., 2019), is employed in our experiments. FESOM2.0 solves the primitive equations in the Boussinesq and hydrostatic approximations. It adopts an unstructured triangular-mesh framework, with scalar degrees of freedom located at vertices and horizontal velocities located at triangle centers. Additionally, the Finite-Element Sea Ice Model (FESIM; Danilov et al., 2015) is incorporated into FESOM2.0 as a set of subroutines. FESIM solves the modified elastic–viscous–plastic (mEVP) dynamical equations, enabling a reduction in subcycling steps while maintaining numerical stability (Kimmritz et al., 2017; Koldunov et al., 2019).

Figure 2 presents the horizontal resolution of the PD and LGM mesh configurations used in this study. The PD mesh consists of 126 858 nodes and 244 659 cells, while the LGM mesh comprises 104 425 nodes and 203 142 cells. Both meshes have the same nominal resolution of 1° for most parts of the global ocean, which translates to a resolution of approximately 25 km north of 50° N, a resolution of about 0.33° at the Equator, and a 10 km resolution for the Arctic Ocean and Bering Sea.

Figure 2 Horizontal resolutions of the (a) PD and (b) LGM mesh configurations used in this study.

The K-profile parameterization (Large et al., 1994) is utilized universally for targeting surface ocean mixing, whereas a constant vertical background diffusivity (k_bg) is employed to manage the effects of various background mixing mechanisms. In FESOM2.0, the default value for k_bg is set to $\mathrm{0.1}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹. Furthermore, the tidal-mixing parameterization by Schmittner and Egbert (2014), drawing on the foundational work of Jayne and St. Laurent (2001) and Simmons et al. (2004), is incorporated. This parameterization uniquely accounts for the influence of subgrid-scale bathymetry on the penetration depth of energy inputs and differentiates between diurnal and semidiurnal tidal effects. The tidal diapycnal diffusivity (k_{v_tidal}) is given by

$$\begin{array}{}\text{(6)}& k_{\mathrm{v}\mathrm{\_}\mathrm{tidal}}={\displaystyle \frac{\mathrm{\Gamma }\mathit{ϵ}}{N^{\mathrm{2}}}},\end{array}$$

where Γ is the mixing efficiency, which is set to 0.2, and N² is the buoyancy frequency. The rate of tidal-energy dissipation (ϵ) is given as

$$\begin{array}{}\text{(7)}& \mathit{ϵ}={\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\sum }_{z^{\prime }>z}^H{\sum }^{\mathrm{TC}}{q}_{\mathrm{TC}}{D}_{\mathrm{IT},\mathrm{TC}}(x,y)F(z,z^{\prime }),\end{array}$$

where $D_{\mathrm{IT},\mathrm{TC}}(x,y)$ is the internal-tide energy flux from barotropic tides to internal tides from the tidal model and F is the vertical decay function with an e-folding depth of 500 m above the seafloor (H). The local dissipation efficiency (q_TC) accounts for the critical latitude (y_c) of diurnal and semidiurnal tidal constituents (TCs):

$$\begin{array}{}\text{(8)}& q_{\mathrm{TC}}=\left\{\begin{array}{l}\mathrm{1},\mathrm{when}\left|y\right|>y_{\mathrm{c},\mathrm{TC}},\\ \mathrm{0.33},\mathrm{otherwise},\end{array}\right.,\end{array}$$

where y_c = 30° for the diurnal constituents (K1 and O1) and y_c = 72° for the semidiurnal constituents (M2 and S2).

3 Model and experiments

In the Methods section, we extract the internal-tide dissipation, denoted as D_IT, from the global tidal model. In comparison to PD values, D_IT for the four principal tidal constituents (M2, S2, K1, and O1) during the LGM shows almost a 3-fold increase, escalating from 1.31 to 3.41 TW (Table 1). The predominant contributor to this shift is the M2 tide, with a period of 12.42 h, which closely matches the North Atlantic basin's period of 12.66 h (Muller, 2008), creating resonance. During the LGM, the removal of continental shelves decreased damping, causing a significant increase in the M2 tide. Its value surged from a PD level of 0.89 TW and reached 2.94 TW during the LGM. These values are in close agreement with previous research findings (Arbic et al., 2004a; Egbert et al., 2004; Green, 2010; Griffiths and Peltier, 2008, 2009; Wilmes and Green, 2014). The horizontal distributions of D_IT (Fig. 1) are used as input for a tidal-mixing parameterization in FESOM2.0.

The experiments are designed to explore how tidal mixing impacts the glacial AMOC, with specifics outlined in Table 4. PD simulations are forced using the 1958–2020 period from the reanalysis dataset (JRA55-do v1.4.0) by Kobayashi et al. (2015), which represents the second global atmospheric reanalysis conducted by the Japan Meteorological Agency (JMA). As of 2009, it has employed the TL319 version of JMA's operational data assimilation system since 1958. The Japanese 55-year Reanalysis (JRA-55) addresses several issues identified in previous reanalyses and enhances the temporal consistency of temperature analysis, making it suitable for examining multidecadal variability and conducting climate simulations. We repeated each PD case simulation using the 1958–2020 JRA-55 data five times to achieve simulation stability. No significant trend was detected, and we utilized the average results from the final cycle for our analysis. The simulations conducted for the PD scenario using FESOM2.0 have been thoroughly validated. Detailed assessments and descriptions of the PD case's configuration can be found in Scholz et al. (2019) and (2022). Figure 3 presents a comparison of Atlantic Ocean temperature and salinity data from our PD case with World Ocean Atlas (WOA) 2018 data. The results indicate that our model accurately reproduces the temperature and salinity structures of the modern ocean. This provides a solid foundation for further simulations of the LGM ocean and the study of the role of tides.

Figure 3 Comparison of salinity and temperature between the WOA 2018 data and the PD simulation with respect to the Atlantic Ocean.

Regarding the LGM simulations, differences in model configuration are attributed to surface forcing and initial conditions. Both of these are taken from the LGMW case in Zhang et al. (2013). It is worth noting that selecting an appropriate forcing for the LGM simulations is crucial. In Knorr et al. (2021), a comparison of the LGM simulations of Atlantic Ocean temperature and salinity structure was conducted among different models from the Paleoclimate Modelling Intercomparison Project (PMIP; Braconnot et al., 2007; Weber et al., 2007) and the Coupled Model Intercomparison Project (CMIP; Braconnot et al., 2012; Kageyama et al., 2017). From these, we selected the simulation results from the well-performing climate model CoSMoS (Zhang et al., 2013) as the surface forcing for the LGM simulations in this study. This provides a solid foundation for accurately simulating the glacial ocean. The LGM simulations are executed over a duration of 600 years to achieve a quasi-equilibrium state. A time series depicting the strength of the AMOC in the LGM cases is presented in Fig. S1 in the Supplement. The concluding 62 years of this period were selected. The simulations are summarized in Table 4.

Table 4 Experimental design of the simulations in this study.

4 Results

The AMOC strength varies between 11.6 and 15.5 Sv (sverdrups) for the PD and between 13.3 and 13.9 Sv for the LGM (Fig. 4; Table S1 in the Supplement). It is noteworthy that both LGM simulations (i.e., the LGM and LGM_tidal simulations) exhibit a shoaled AMOC, identified at approximately 1700 m depth. This suggests that the inclusion of tidal mixing does not affect the glacial AMOC's configuration. Instead, stratification is identified as a crucial factor contributing to the shoaled AMOC, which is more pronounced in the upper and middle regions of the glacial ocean. Stratification, on the one hand, signifies the buoyancy forces encountered by water masses during their descent. On the other hand, it exerts a notable influence on the model's tidal diffusivities, as outlined in Eq. (6). Our tidal-model results indicate that internal-tide dissipation (D_IT) in the North Atlantic during the LGM reached 0.81 TW, a significant increase compared to the value of 0.13 TW observed in the PD, representing a 6-fold increase. Additionally, during the LGM, the average squared buoyancy frequency (N²) at depths of 1000–3000 m in the Atlantic corresponds to $\mathrm{1.89}\times {\mathrm{10}}^{-\mathrm{6}}$ s⁻², compared to $\mathrm{1.02}\times {\mathrm{10}}^{-\mathrm{6}}$ s⁻² in the PD. These changes do not suggest tidal dominance but rather imply a stronger vertical tidal diffusivity (shown in Fig. 4b and e). This enhanced vertical mixing aligns with the observed reductions in the vertical gradients of radiocarbon and δ¹³C in the deep Atlantic during the LGM (Skinner et al., 2017; Muglia et al., 2018; Peterson et al., 2014; Molina-Kescher et al., 2016; Sikes et al., 2016).

Figure 4 Shown are the AMOC (left) as well as zonally averaged distributions for the squared buoyancy frequency (middle) and tidal diffusivity (right) across the Atlantic Ocean. Simulations are as listed in Table 4.

However, the incorporation of tidal-mixing processes results in a substantial increase in the generation of AABW, with a magnitude of −7.9 Sv, significantly exceeding the PD estimate of −4.2 Sv. This enhancement is in agreement with the results derived from paleoproxy data, which indicate intensified glacial AABW (Curry and Oppo, 2005; Zhang et al., 2017). Consequently, while LGM tides may not modify the AMOC, they exert a significant influence on AABW formation. Thus, accounting for tidal effects remains essential in conducting climate modeling for the LGM period.

In the PD_tidal experiment, incorporating the tidal-mixing parameterization does not alter the geometry of the AMOC. However, the PD_glacial_tidal simulation, which includes LGM tidal dissipation, reveals distinct dynamics. This simulation demonstrates a significant increase in both the strength and depth of the AMOC, extending to near-benthic layers at approximately 35° N (Fig. 4c and Table S1), and is accompanied by notably reduced stratification. These observations suggest that the effects and dynamics of enhanced tidal dissipation differ substantially under varying ocean stratification intensities during the LGM and PD periods.

To investigate the origins of various stratifications from the perspectives of temperature and salinity, we conduct a further analysis of the temperature and salinity distributions in the Atlantic Ocean across different cases (Fig. 5). Under PD forcing, the surface salinity below a latitude of 40° N is significantly higher than that in the middle and lower layers, which weakens ocean stratification. This sharply contrasts the high abyssal salinity observed during the LGM, a defining feature of the glacial ocean (Adkins et al., 2002; Knorr et al., 2021). Regarding temperature, a decrease from the surface to the seabed is observed in both the PD and LGM scenarios. Notably, during the LGM, the simulated temperature exhibits a pronounced decrease above a depth of 2 km, decreasing to 0 °C at this level. Below a depth of 2 km, the ocean is relatively homogeneous, with temperatures close to the freezing point, indicating a cold and well-mixed deep ocean, consistent with paleoproxy data (Adkins et al., 2002). This steeper temperature gradient enhances stratification in the upper layers of the glacial ocean. Collectively, the high abyssal salinity and swift vertical temperature decline significantly contribute to more pronounced stratification within the glacial ocean. Accurately replicating these temperature and salinity features is crucial in the climate modeling of the LGM.

Figure 5 Zonally averaged distributions across the Atlantic Ocean for the simulations listed in Table 4. The distributions for the squared buoyancy frequency (left), salinity (middle), and potential temperature (right) are shown.

In our comparable consideration of enhanced tidal mixing during the LGM, the discrepancy between our results and those of previous studies (Schmittner et al., 2015; Wilmes et al., 2019) is attributed to the fact that these studies do not reproduce high abyssal salinity and increased stratification in the LGM Atlantic. These features counteract the impact of stronger tidal mixing on the AMOC.

5 Discussion

Tides play a pivotal role in climate dynamics, e.g., by facilitating the release of iceberg armadas during Heinrich events (Arbic et al., 2004b) and serving as the primary driving force behind both vertical and horizontal ice sheet movements at the corresponding marine peripheries (Padman et al., 2018). Focusing on the LGM period reveals that glacial tidal dissipation was approximately 3 times greater than present levels. This increase, combined with the closure of the Bering Strait (Hu et al., 2010), led to reduced freshwater transport to the Atlantic, ostensibly resulting in a strengthened glacial AMOC. However, paleoclimatic proxy data indicate a significant shoaling of the AMOC during the LGM, with an estimated reduction in depth of about 1000 m compared to contemporary conditions (Burke et al., 2015; Lund et al., 2011). The primary aim of this study is to identify the reasons for the shoaled glacial AMOC given the complex interplay of these factors.

Our results indicate that the integration of additional tidal-mixing parameterizations does not significantly influence the AMOC in either the PD or the LGM scenarios. In the PD scenario, the relatively weak tides can be adequately accounted for by background diffusivity (k_bg), thereby negating the necessity for an additional tidal parameterization. Furthermore, during the LGM, tides are unlikely to have played a substantial role in influencing the glacial AMOC due to pronounced ocean stratification. On the one hand, this stratification hampers the mixing of water masses, while on the other hand, it leads to a decrease in the effectiveness of tidal mixing. This is because the buoyancy frequency, which appears in the denominator of the tidal-mixing parameterization (as detailed in Eq. 6 in the Method section), suggests that stronger stratification significantly reduces the impact of tidal dissipation. However, in the abyssal ocean, where there is relatively weak stratification, the pronounced tidal dissipation during the LGM notably enhances the formation of AABW.

The tidal-mixing parameterization only considers locally dissipated energy, which only accounts for one-third of the total energy (Jayne, 2009; Schmittner and Egbert, 2014). The remaining two-thirds of the energy is dissipated in the far-field, where background diffusivity (k_bg) is employed to represent this dissipation. Consequently, we calculated the far-field dissipation due to k_bg and the local tidal dissipation due to k_{v_tidal} for each simulation using the Osborn (1980) formula, expressed as follows:

$$P=\int \mathit{\rho }\mathit{ϵ}\mathrm{d}V={\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }}}\int \mathit{\rho }k{N}^{\mathrm{2}}\mathrm{d}V.$$

The results are presented in Table 5. Additionally, we conducted another experiment, LGM_tidal_3, in which k_bg increased from $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ to $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹ for comparison. The results indicate that the LGM_tidal experiment underestimates tidal energy, reaching only 2.05 TW. In contrast, the LGM_tidal_3 experiment shows no such underestimation, with its energy reaching 3.92 TW.

Table 5 Summary of the energy consumption due to diapycnal mixing.

Figure 6 presents the AMOC geometry for the LGM_tidal and LGM_tidal_3 experiments. The geometry of the AMOC in LGM_tidal_3 experiment remains relatively shallow, without significant changes, which further supports our study's conclusions. The only notable change is in the AMOC strength, which increased from 13.3 to 15.4 Sv. This underscores the necessity of employing the tidal-mixing parameterization and the importance of appropriately adjusting the background diffusivity (k_bg). Additionally, this result is the first to demonstrate that the shallower geometry of the glacial AMOC during the LGM remains unchanged, even when accounting for far-field tidal dissipation.

Figure 6 AMOC stream functions (measured in sverdrups) for the LGM_tidal and LGM_tidal_3 simulations.

Applying enhanced LGM tidal dissipation to PD conditions (in the case of the PD_glacial_tidal simulation), where ocean stratification is significantly weaker than that under LGM conditions, yields entirely different results. In this scenario, the amplified tidal dissipation induces a deeper and more potent AMOC. We propose a potential positive feedback mechanism that accounts for the increased AMOC during termination. A reduced ocean stratification during the initial phase of termination enhances the effectiveness of tidal mixing, a process analogous to the one discussed above. This increased tidal mixing will affect the ocean more efficiently, further weakening ocean stratification, thereby increasing tidal mixing. This initiates a positive feedback loop that culminates in reduced stratification and a more vigorous and deeper AMOC.

6 Conclusions

The concept of enhanced glacial-ocean stratification, which potentially results from the cooling and salinification of glacial AABW and could lead to a shoaled AMOC during the LGM, has previously been discussed (Jansen and Nadeau, 2016; Jansen, 2017; Klockmann et al., 2016). However, until now, no research has directly demonstrated a shoaled AMOC under real LGM forcing conditions, including the impact of increased glacial tidal dissipation. Montenegro et al. (2007) proposed that LGM tides have a minimal impact on the AMOC, attributing this to a potential underestimation of tidal dissipation during the LGM. In contrast, Schmittner et al. (2015) and Wilmes et al. (2019) suggested a significant enhancement and deepening of the North Atlantic overturning cell under vigorous glacial tidal dissipation. It is noteworthy that Wilmes et al. (2021) obtained a relatively shoaled LGM AMOC through the artificial reduction of meridional moisture flux and precipitation at high latitudes.

Our study is the first to directly demonstrate a stratified ocean and a shoaled AMOC under real LGM forcing conditions, without any artificial modifications, despite the presence of increased glacial tidal dissipation. We suggest that accurately simulating the high salinity of the deep sea and the rapid temperature changes in the ocean's upper layers is crucial for correctly reproducing a stratified glacial ocean. In such an environment, the significant tidal dissipation during the LGM was not sufficient enough to counter the increased ocean stratification, leading to a shoaled AMOC. Furthermore, we emphasize that this notable glacial tidal dissipation plays a critical role in strengthening AABW during the LGM.

Our results highlight the dominance of background conditions and mixing with respect to ocean circulation dynamics, as well as possible complex feedbacks in the Earth system (Lohmann et al., 2020). Here, we use an ocean-only model, which means the LGM atmospheric forcing is kept constant. Consequently, this approach has a limitation: it does not account for interactions between the ocean (or sea ice) and the atmosphere. As a logical next step, we plan to incorporate tidal-energy dissipation into fully coupled Earth system models (Zhang et al., 2014; Liu et al., 2009) to further explore AMOC dynamics during deglaciation.

Appendix A: Detailed iterative process for eliminating N² sensitivities between the tidal model and FESOM2.0

To address the N² sensitivities and interactions between FESOM2.0 and the tidal model, we employed an iterative process for the PD_tidal and LGM_tidal simulations. Here, we provide a detailed description of this process, employing the LGM_tidal simulation as an example.

The steps of the iterative process are as follows:

1. For the initial input, we begin by obtaining N² from the LGM simulation, which does not include tidal mixing. This initial N² is then used as input for the tidal model.

2. In the first iteration (LGM_tidal1), the tidal model, which uses the initial N², calculates the tidal dissipation, which is then fed back into FESOM2.0, producing the first experimental result (LGM_tidal1).

3. In the second iteration (LGM_tidal2), N² from the LGM_tidal1 simulation is inputted into the tidal model once again to generate a new tidal dissipation, which is incorporated into FESOM2.0, and the model is run to obtain LGM_tidal2.

4. The final output (LGM_tidal2) represents the second iteration and is the LGM_tidal experiment presented in our paper.

Figure A1 Changes in depth-averaged vertical N² values across iterations for the LGM and PD simulations.

Fig. A1 illustrates the changes in depth-averaged vertical N² values during these iterations for both the PD and LGM scenarios. For the LGM simulations, the primary change from the initial LGM iteration to the LGM_tidal1 iteration involved a decrease in N² in the Arctic, and from the LGM_tidal1 iteration to the LGM_tidal2 iteration, minimal changes were observed. However, for the PD simulations, no significant changes in N² were observed throughout the iterations, effectively minimizing the mutual influences between N² in the tidal model and that in the OGCM after the iterative process. This iterative approach ensures the stability and accuracy of our model results, reducing the sensitivity of N² to tidal-dissipation feedbacks.