Late Pleistocene Evolution of Tides and Tidal Dissipation

Studies of the Last Glacial Maximum (LGM; 26.5–19 ka) tides showed strong enhancements in open ocean tidal amplitudes and dissipation rates; however, changes prior to the LGM remain largely unexplored. Using two different ice sheet and sea level reconstructions, we explicitly simulate the evolution of the leading semi‐diurnal and diurnal tidal constituents (M2, S2, K1, and O1) over the last glacial cycle with a global tide model. Both sets of simulations show that global changes, dominated by the Atlantic, take place for the semi‐diurnal constituents, while changes for the diurnal constituents are mainly regional. Irrespective of the reconstruction, open ocean dissipation peaks during the sea level lowstands of MIS 2 (∼20 ka) and MIS 4 (∼60 ka), although dissipation values prior to MIS 2 are sensitive to differences in reconstructed ice sheet extent. Using the statistically significant relationship between global mean sea level and dissipation, we apply regression analysis to infer open ocean and shelf dissipation, respectively, over the last four glacial cycles back to 430 ka. Our analysis shows that open ocean tidal energy was probably increased for most of this period, peaking during glacial maxima, and returning to near‐present‐day values during interglacials. Due to tidal resonance during glacial phases, small changes in bathymetry could have caused large changes in tidal amplitudes and dissipation, emphasizing the need for accurate ice margin reconstructions. During glacial phases, once global mean sea level decreased by more than ∼100 m, the amount of open ocean tidal energy available for ocean mixing approximately doubled.


Introduction
Tides are important for numerous processes in the ocean: in coastal areas, they shape intertidal ecosystems and morphology.In shelf-sea areas, tidal dynamics determine the location of tidal mixing fronts, which separate seasonally stratified waters from year-round mixed waters (e.g., Simpson & Pingree, 1978).This partitioning is important for shelf sea ecosystems, biogeochemical cycles and the export of CO 2 from the shelf seas to the deep ocean (Thomas et al., 2004).In the open ocean, tidally driven mixing supplies approximately half (∼1 TW) of the energy necessary to sustain the large-scale meridional overturning circulation (Ferrari & Wunsch, 2009;Wunsch & Ferrari, 2004).Recent work (Wilmes et al., 2021) demonstrated that, for the Last Glacial Maximum (26.5-19 thousand years before present (ka); the LGM hereafter), strong increases in tidal mixing (due to the sea level lowstand and associated changes in tidal dynamics) could be constrained from sediment carbon isotopes.As sea level index points (SLIPs) in coastal areas are generally related to a given high or low tide level and not mean sea level, knowledge of past tidal range changes is important for reconstructing past sea levels (see e.g., Ward et al., 2016).Furthermore, marine terminating ice sheet dynamics are affected by tidal dynamics which influence grounding line movement (Batchelor et al., 2023;Milillo et al., 2017), basal melting (Anselin et al., 2023;Milillo et al., 2017), ice shelf flexure (Walker et al., 2013) and ice ow (Anandakrishnan, 2003;Bindschadler et al., 2003;Gudmundsson, 2007).
Reconstructions of global tides and tidal dissipation during the Quaternary (2.58 Ma to present) have generally focused on the last ∼25 thousand years (kyr) encompassing the LGM, the deglacial (19-11.7 ka), and the Holocene (11.7 ka to present) (Egbert et al., 2004;Green, 2010;Griths & Peltier, 2008, 2009;Sulzbach et al., 2023;Uehara et al., 2006;Wilmes et al., 2019Wilmes et al., , 2022;;Wilmes & Green, 2014).These investigations showed surprising results: the tides were strongly enhanced in the Atlantic during the LGM, especially in the semi-diurnal band, with tidal energy dissipation (i.e., the loss of energy of the tide to bed friction and to the internal tide) for the M 2 tide a factor 2-3 larger than at present in the open ocean.Changes in the North Atlantic are thought to have been particularly strong with amplitudes tripling with respect to present and exceeding 6 m in the Labrador Sea during the LGM (Griths & Peltier, 2008, 2009;Wilmes & Green, 2014).These amplifications resulted from changes in ocean basin shape driven by the ∼130 m global mean sea level (GMSL) drop and associated increases in ice sheet extent.Together, these factors rendered the Atlantic more conducive to resonant amplification of the semi-diurnal tides, thus leading to substantial increases in amplitudes and tidal energy dissipation (e.g., Egbert et al., 2004;Green, 2010).Through the deglacial and the Holocene, energy losses in the deep ocean decreased dramatically whilst the shelf seas (which were emerged during the LGM) re-flooded and became more tidally energetic.
On longer time scales, covering the middle and late Pleistocene (∼770-11.7 ka), sea level fluctuated by 130-145 m as climate cycled between glacial and interglacial phases (e.g., Fox-Kemper et al., 2021).During this period, ice sheet extent followed a saw-tooth shaped pattern: glacial phases were generally characterized by a long-term gradual cooling culminating in a glacial maximum with peak ice sheet extent and sea level lowstands.Subsequently, climate warmed rapidly and transitioned to interglacial conditions, with climate similar to, or warmer than the pre-industrial.GMSL attained highstands of up to +15 m and ice sheet extent was similar to, or less than present.The last glacial period spanned from the end of the Last Interglacial (LIG; 130-115 ka; Marine Isotope Stage (MIS) 5e) to the onset of the Holocene, with significant variations in global mean sea level (Lambeck, 2004;Siddall et al., 2003;Spratt & Lisiecki, 2016;Waelbroeck et al., 2002) (see Figure 1a).During the early part of the glacial period , GMSL fluctuated between relative high-and lowstands of −9 and −50 m, respectively (e.g., Creveling et al., 2017).During MIS 4 (71-57 ka), a GMSL lowstand of around −80 m was reached.Thereafter, GMSL rose to a relative highstand during MIS 3 (57-29 ka), though the exact magnitude remains debated (e.g., Dalton et al., 2022).The GMSL lowstand of ∼−130 m during the LGM (MIS 2) was reached between 26.5 and 20 ka (Clark et al., 2009;Gowan et al., 2021;Lambeck et al., 2014;Peltier et al., 2015).At the onset of the deglacial period (19-11.7 ka), GMSL first rose gradually, then more rapidly by around 100 m until the early Holocene (∼8 ka) when present-day levels were reached.Because tides are sensitive to water depth and basin geometry, it is expected that they are affected by these large sea level changes, as previously seen for the period spanning from the LGM to the present (e.g., Egbert, 2004;Green, 2010;Sulzbach et al., 2023;Uehara et al., 2006;Wilmes et al., 2022;Wilmes & Green, 2014).
In this work, we aim to extend our knowledge of Pleistocene tidal dynamics back to ∼430 ka.First, we explicitly model tides over the past glacial cycle covering the period from the LIG to present using two different ice sheet and sea level reconstructions, thus expanding the work in Wilmes et al. (2022) to the entire last glacial cycle.Second, because no spatially and temporally highly resolved global sea level and ice sheet reconstructions exist for the previous multiple glacial cycles, we use the last glacial cycle simulations together with uniform sea level change simulations and extrapolate tidal energy dissipation back to 430 ka based on linear regression analysis.We aim to (a) produce global spatially-varying dissipation estimates for the last glacial cycle which can be used for modeling of the late Pleistocene climate and ocean circulation; and (b) improve our understanding of late Pleistocene tidal dynamics which is relevant for, for example, SLIPs, ice sheet dynamics, or shelf sea oceanographic processes.
where U is the depth integrated volume transport, H denotes water depth, f is the Coriolis vector, g is the gravitational constant, ζ denotes tidal elevation, ζ EQ stands for the equilibrium tidal elevation, and ζ SAL is the tidal elevation due to self-attraction and loading (SAL; i.e., the combined effects of gravitational attraction among the water masses, seafloor deformation, and associated changes in the gravity potential).F = F B + F IT represents frictional losses due to bed friction (F B ) and tidal conversion (F IT ).The former is represented by the standard quadratic law: where C d = 0.003 is a drag coefficient, and u tidal velocity.The energy losses to the internal tide, F IT = C IT U, depend on a conversion coefficient C IT given by (Green & Huber, 2013;Zaron & Egbert, 2006) where γ = 37.5 is a scaling factor (see Zaron & Egbert, 2006, for more details), N b is the buoyancy frequency at the sea-bed,  N is the vertical average of the buoyancy frequency, and ω is the frequency of the tidal constituent under evaluation.γ was tuned following the process described in Wilmes and Green (2014) to both minimize present-day amplitude root-mean square errors against TPXO9 and t TPXO9 dissipation values.We assume horizontally uniform abyssal stratification by parameterizing the buoyancy frequency N through N(z) = N 0 e (−z/1300) with N 0 = 5.24 × 10 −3 .Explorations in Wilmes and Green (2014) and Schmittner et al. (2015) showed that tidal dissipation is rather insensitive to glacial-interglacial stratification changes.Throughout the paper, we calculate tidal dissipation rates associated with the combined action of F B and F IT via the energy balance method outlined in Ray et al. (2003).

Last Glacial Cycle
In our tide model runs, we represent the ice sheet and sea level history over the last glacial cycle with two different ice sheet reconstructions, one external reconstruction by Gowan et al. (2021) and one derived by us, using a global gravitationally self-consistent sea level model.For the latter, we simulate sea level based on the ICE6G_C ice history (Argus et al., 2014;Peltier et al., 2015) that spans the period 122 ka to present with a temporal resolution of 250-500 years prior to 21 ka, 1,000 years up to 32 ka, and 2,000 years thereafter (referred to as ICE-6G from this point onwards).In the sea level model, we compute the time-varying deformation of a rotating, Maxwell viscoelastic Earth model with a depth-dependent Earth structure (e.g., Kendall et al., 2005;Milne & Mitrovica, 1996).For consistency with the global ICE6G_C ice model, we adopt the corresponding VM5a depth-dependent Earth model to represent global Earth structure (Argus et al., 2014;Peltier et al., 2015).Importantly, for our application here, the sea level model includes migrating shorelines and the inundation of water into regions previously covered by marine ice (Mitrovica & Milne, 2003).The model formulation is solved using a pseudo-spectral numerical scheme up to spherical harmonic degrees and order 256 (Kendall et al., 2005).
As an alternative to the boundary conditions derived from the ICE-6G_C ice sheet history, we use the ice sheet and sea level reconstruction from Gowan et al. (2021).Gowan et al. (2021) presented a global high-resolution ice sheet reconstruction for the period 80 ka to present.The reconstruction is consistent with ice physics but was inferred independently of far-field sea level and δ 18 O proxy records (i.e., based on ice sheet margins and constructed to be consistent with simple ice sheet physics).The sea level predictions provided by Gowan et al. (2021) also rely on a global gravitationally self-consistent sea level model (SELEN; e.g., Spada & Stocchi, 2007) with a depth-dependent Earth structure, as well as shoreline migration.These reconstructions have a lower temporal resolution of 2,500 years.
For our study, both ice sheet reconstructions and the associated sea level change fields were interpolated to the finite-difference grid of OTIS (1/8° spacing in both latitude and longitude) and added to the present-day base topography RTopo-2 (Schaffer et al., 2016).We note that both the ice sheet history and the solid Earth structure, and thus sea level, differ between the two approaches.For the ICE-6G sea levels, tide runs were performed at 2,000-year intervals.Each simulation was run with M 2 , S 2 , K 1 and O 1 equilibrium tidal forcing and a simplified SAL scheme that sets ζ SAL = βζ, with β = 0.1.The computational grid extends all the way to 89°N, where it is bounded by an artificial vertical wall eliminating the need for open boundaries.These two simplifications (scalar SAL and pole cap) potentially suppress resonant behavior of the Arctic Ocean during the LGM (Griths & Peltier, 2009;Sulzbach et al., 2023).However, the Arctic Ocean is of secondary importance for globally integrated dissipation rates, and simulations with a more accurate iterative SAL scheme (Wilmes et al., 2022) lend credence to our choice of β = 0.1, at least for the deglacial.
Additionally, to separate the effect of GMSL changes from land-ocean boundary changes and effects of non-uniform sea level change, we conducted tidal simulations with uniform sea level changes.Here, sea level was uniformly changed from −135 m to +20 m in steps of 5 m.Again, RTopo-2 at 1/8° × 1/8° grid spacing was used as the present-day base topography and the simulations were carried out for M 2 and K 1 .
The present-day control simulation (i.e., for 0 ka) was benchmarked against the global tidal solution TPXO9-atlas-v5 (https://www.tpxo.net/global/tpxo9-atlas)by calculating amplitude root-mean square errors and comparing globally integrated dissipation values; see Table 1.Evidently, the model yields realistic solutions for both diurnal and semi-diurnal tidal constituents.

Multiple Glacial Cycles (430 ka-present)
Owing to the lack of global bathymetry, ice sheet and sea level reconstructions spanning multiple glacial cycles, it is currently not possible to explicitly model tidal dynamics prior to the LIG.We therefore take a different, novel approach to infer tidal dissipation prior to the last glacial cycle to gain a first-order understanding of tidal dynamics during this period.We first establish regression models between GMSL and globally integrated deep ocean and shelf sea dissipation rates, respectively, using the three simulation sets described above, that is, runs for (a) the last glacial cycle based on ICE-6G, (c) the last glacial cycle based on Gowan et al. (2021), and (c) uniform GMSL changes.Using the global sea level reconstruction by Spratt and Lisiecki (2016) (thereafter SL16), each regression model is then used to infer tidal dissipation for open ocean and shelf ocean dissipation back to 430 ka.Because GMSL is predominantly driven by global ice volume changes, and the ice sheet volume changes over the last multiple glacial cycles generally follow similar patterns (e.g., waxing and waning of the Laurentide and Fennoscandinavian ice sheets) (e.g., Batchelor et al., 2019), we here use GMSL as a proxy for the combined spatially-varying ice volume and sea level changes.The different bathymetry reconstructions provide a measure of the uncertainty introduced by the differing ice margins.
For the regression models, we chose a polynomial regression of order two between GMSL (sl) and globally integrated tidal dissipation (D) for the deep and shelf ocean (D deep and D shelf ), respectively, for each tidal constituent.The quadratic t accounts for non-linear interactions between GMSL and dissipation due to, for example, resonance effects.The relationship between GMSL and tidal dissipation thus takes the following form.
where {β d0 , β d1 , β d2 } and {β s0 , β s1 , β s2 } represent the regression coefficients for deep and shelf dissipation for each set of model simulations, respectively.We calculate regression coeficients for the relationships between GMSL and deep and shelf dissipation for all three sets of tide runs, respectively.If the relationship between sl 2 and D is not significant at the 95% confidence level, the order of the polynomial is reduced to one (i.e., linear regression).
The relationship between sl and D deep and D shelf , respectively, is calculated separately because, in general, deep and shelf dissipation behaved in an anti-correlated manner in relation to GMSL (i.e., when GMSL decreases, open ocean dissipation increases and shelf dissipation decreases).The globally integrated dissipation D total is given by the sum of D deep and D shelf .
Using each regression model and the SL16 GMSL reconstruction, we then calculate time series of late Pleistocene deep, shelf, and total dissipation.Dissipation values are predicted for each constituent (M 2 , S 2 , K 1 , and O 1 ) as per Equations ( 5)-( 7) and associated standard deviations (σ) are deduced by applying the laws of variance propagation.
where sl denotes GMSL from the SL16 reconstruction with the provided formal error σ sl , and σ β… denotes the standard deviations associated with the respective regression coefficients.

Ice Sheet and Sea Level Evolution
The two ice sheet and sea level reconstructions show pronounced differences during large parts of the last glacial cycle (Figures 1a and 2).In the following, we will focus on the differences between the two.For a detailed  Ice Sheet is slightly more extensive than in Gowan, with the Laurentide Ice Sheet and the Cordilleran Ice Sheet joined up, whereas they remain separated in Gowan.However, the Gowan reconstruction has more extensive ice around Antarctica.Overall, as a result of these differences, ICE-6G sea levels show a slightly lower and earlier MIS 4 GMSL lowstand than Gowan.However, whilst GMSL leading into MIS four is similar, there are large regional differences in sea level due to the differences in ice history.The largest differences in ice sheet extent and sea level between the two reconstructions are seen during MIS 3. ICE-6G shows extensive Northern Hemisphere ice sheets throughout MIS 3 whereas in Gowan, ice sheet extent is strongly reduced around 40 ka, with the Fennoscandinavian and Cordilleran Ice Sheets mostly melted and the Laurentide Ice Sheet strongly reduced.These discrepancies lead to an offset in GMSL of around 60 m in the middle of MIS 3 and shelf sea area is approximately doubled in Gowan in comparison to ICE-6G.Toward the LGM, both reconstructions show expanding ice sheets and a drop in GMSL.ICE-6G suggests slightly more extensive Northern Hemisphere ice sheets, whereas Gowan shows a stronger expansion of the Antarctic ice sheet margins.The ICE-6G sea level lowstand occurs around 26 ka, whereas, for Gowan, it does not take place until 20 ka, and is around 20 m less than for ICE-6G.From 20 ka to present, both simulations show a similar GMSL evolution.However, again, there are pronounced differences in regional sea levels (locally, >50 m offsets) due to differences in instantaneous ice sheet loading but also ice history.

Amplitude Evolution
Large changes in tidal amplitudes (regionally in excess of 4 m) occur over the last glacial cycle for both sea level reconstructions, but, owing to the discrepancies in sea level and ice sheet history, there are pronounced offsets between the tidal histories (see Figures 3 and 4 for total amplitude changes).In the following, we discuss how each constituent contributes to the total amplitude changes, focusing the description on the areas that show the largest changes and laying the emphasis on M 2 , which shows the largest absolute changes.amplitudes increase further to peak round 20 ka, however, Weddell Sea amplitudes are much reduced in comparison to MIS 3 because of the more extensive grounded ice in the bay.During the deglacial, pronounced drops in North Atlantic and Arctic M 2 amplitudes take place, and by 8 ka tidal amplitude reach near-present day levels.
For the Gowan simulations, M 2 amplitudes (Figure S5 in Supporting Information S1) between 80 ka and 72.5 ka show similar levels to present, apart from considerably larger amplitudes in the Weddell and Ross Sea area.As the sea level lowstand around 60 ka is approached, North Atlantic amplitudes strongly increase.These enhancements persist (albeit at a smaller magnitude) until ∼55 ka, after which the North Atlantic oscillates between periods of larger and smaller M 2 amplitudes with peaks around 45 ka, 37.5 ka and 32.5 ka.Notable are also strongly increased tides in the Ross Sea between 50 and 40 ka.Large North Atlantic amplitudes develop toward the LGM and during the early to mid deglacial (peaking around 17.5 ka and persisting until 12.5 ka).During LGM, strong enhancements in Arctic tides can also be seen.
For the ICE-6G runs, S 2 amplitudes (Figure S2 in Supporting Information S1) are similar to present at the end of the LIG.During the early glacial, Labrador Sea and North Atlantic amplitudes increase, but the increases are less pronounced than for M 2 .Furthermore, large S 2 amplitudes occur in the Coral Sea (NE Australia).The two sea level lowstands (∼60 ka and LGM) are characterized by large North Atlantic tides and strong enhancements in Arctic S 2 amplitudes.S 2 amplitudes in the Gowan simulations (Figure S6 in Supporting Information S1) follow a similar picture as M 2 amplitudes albeit with a lower magnitude.Notable is that Arctic S 2 amplitudes enhancements are much reduced in the Gowan simulations and are only present during the middle of the LGM.
In the ICE-6G runs, K 1 and O 1 tidal changes (Figure S4 in Supporting Information S1) are mainly regional.
During the sea level lowstands around 62 ka and during the LGM, Pacific shelf seas (Sea of Okhotsk, South China Sea and Banda Sea) become resonant.During the remainder of the last glacial cycle, amplitudes remain close to their present-day levels.
Similarly, both K 1 and O 1 amplitudes are characterized by mainly regional changes in the Gowan simulations (Figures S7 and S8 in Supporting Information S1).During periods of sea level lowstands (∼62.5-55ka and 27.5-15 ka), increased amplitudes can be seen around Antarctica and in the shelf seas of the Pacific, whereas other parts of the ocean seem rather insensitive to the large sea level and ocean basin shape changes.

Dissipation Evolution
For M 2 , global open ocean dissipation in the ICE-6G simulation is at near-present-day values during the LIG (Figure 1b).It approximately doubles over a period of ∼10 kyr between 120 and 110 ka to ∼2 TW when North and South Atlantic dissipation increases (see Figure 5).It remains elevated at this level until around 70 ka when dissipation rises by a further 25% (0.5 TW) driven by Arctic and North and South Atlantic dissipation increases.
Open ocean dissipation remains at these elevated values peaking around 65 ka, 40 ka and during the LGM, with only a small dip occurring around 28 ka.From 68 ka to 30 ka, open ocean dissipation in the Southern Ocean more than doubled in the ICE-6G run, owing to increases in the Weddell Sea area which enhance Southern Hemisphere dissipation to levels greater than in the Northern Hemisphere (Figure S9 in Supporting Information S1).During MIS 4 to MIS 2, dissipation increases around the European Shelf, in the Labrador Sea, Denmark Strait, Norwegian Sea, along the mid-Atlantic ridge, in the Arctic Basin and also in the South Atlantic.In contrast, the Gowan simulations show near present-day open ocean dissipation values prior to 60 ka, with contributions from the Southern Hemisphere (particularly around Antarctica) accounting for the small observed changes (Figure S9 in Supporting Information S1).During sea level lowstand at ∼60 ka, Gowan dissipation values increase by ∼50% through Atlantic enhancements (see Figure 6).Thereafter, dissipation decreases slightly, but oscillates through three distinct troughs and peaks.It is notable that Hudson Bay & Strait dissipation is significantly anti-correlated with North Atlantic dissipation, suggesting that when Hudson Bay & Strait lies dry or is ice covered, North Atlantic dissipation peaks.Toward the LGM, dissipation increases again, and reaches its peak during the early deglacial phase with a 90% increase relative to present driven by increases in Atlantic (especially North Atlantic) and Arctic dissipation, whereas in the ICE-6G simulations, dissipation in the Southern Hemisphere decreases in contrast to the Gowan simulations.
Despite the large differences in shelf area (see Figure 1a), the evolution of shelf sea dissipation is relatively similar between the two sets of simulations.Between 120 and 70 ka, shelf dissipation is at levels similar to present.With the decrease in sea level toward the lowstand around 60 ka, shelf dissipation decreases by ∼40% and remains reduced until the onset of the Holocene with values halved with respect to present during the LGM.In the Gowan simulation, around 40 ka, the values increase to present-day levels and then decrease again to their LGM minimum.Overall, for M 2 , the simulations suggest that open ocean dissipation was enhanced for most of the last glacial period, but that the exact magnitude is dependent on the sea level and ice sheet evolution.Less M 2 energy was lost in the shelf seas between 70 ka and the onset of the Holocene and more M 2 energy dissipated in the open ocean.Overall, the results from the ICE-6G simulations suggest that total dissipation was larger than present for most of the glacial cycle, on average 25% (maximum ∼40%), whereas for the Gowan simulations which only extend back to 80 ka, the average increase was only 7% (maximum 28%).
For S 2 (Figure 1c), in contrast, total dissipation is lower than present for most of the glacial cycle, except for the Holocene (for ICE-6G, the average decrease is 15%).This is driven by pronounced decreases in shelf sea  dissipation which are strongest between 70 ka and the mid-deglacial.The magnitude of the decrease (up to 60% with respect to present) is slightly larger in the ICE-6G simulations due to the larger reduction in shelf area.Deep dissipation is slightly elevated (20%-30%) in the ICE-6G simulation, but reduced with respect to present (apart from the LGM period) in the Gowan simulation.
For K 1 (Figure 1d), the overall pattern is similar to S 2 , with total dissipation being slightly reduced for most of the glacial cycle (for ICE-6G by ∼10%) driven by lower levels of shelf sea dissipation.Again, open ocean dissipation in the ICE-6G simulation is enhanced between 70 ka and the onset of the Holocene and shelf sea dissipation remains reduced; whereas in the Gowan simulation, there are two distinct peaks situated around the sea level lowstands at 60 ka and the LGM.Between the lowstands, open ocean dissipation and shelf sea dissipation returns to present-day values.
For O 1 (Figure 1e), the shelf sea dissipation signal is similar to S 2 and K 1 , but the decreases are compensated by increases in open ocean dissipation, such that total dissipation remains at values near present for most of the glacial cycle and an increase of up to 38% in total dissipation occurs during the peak of the LGM (cf.Sulzbach et al., 2023).
Overall, this means that, for the Gowan simulations, globally integrated dissipation, on average, was close to present-day levels between 80 ka and present, but for the ICE-6G simulations, total dissipation over the whole glacial cycle was on average ∼15% larger than at present (Figure 1f).Deep dissipation in the ICE-6G runs was 57% greater on average (80% larger between 66 ka and 16 ka), and raised by 34% in the Gowan simulation.For shelf dissipation the mean reductions for ICE-6G and Gowan are 22% and 24%, respectively.

Relationship Between GMSL and Tidal Dissipation
For the ICE-6G simulations, open ocean dissipation is significantly anti-correlated with GMSL for all constituents (r < −0.88 for all constituents) that is, when sea level drops, open ocean dissipation increases (see Figure 7).For shelf dissipation, similarly high magnitude significantly positive correlations emerge (r > 0.95) but the relationship is opposite, that is, decreasing sea levels reduce shelf dissipation.In contrast, total dissipation shows a weaker albeit still significant relationship with GMSL, owing to the opposing relationships of open ocean and shelf dissipation.For M 2 and O 1 , total dissipation is anti-correlated with GMSL (i.e., larger dissipation with lower GMSL), whereas for S 2 and K 1 we observe the reverse.For the Gowan simulations, correlations generally have the same sign as the correlations in the ICE-6G simulations, but with slightly reduced magnitudes.Dissipation from the uniform sea level drop simulations shows similarly high negative correlations between GMSL and open ocean dissipation as the ICE-6G simulations, and highly positive correlations between GMSL and shelf dissipation.Next, we evaluate the polynomial regression t between GMSL and both shelf and open ocean dissipation for each tidal constituent for all three sets of simulations, respectively (see Figure 8).For open ocean dissipation, all constituents apart from S 2 show r > 0.95 between actual and estimated dissipation values, and for shelf dissipation all regression models can explain more than 89% of the variability.For the Gowan simulations, the regression t is good for the diurnal constituents but shows a slightly less good t for the semi-diurnal constituents (r = 0.65 and r = 0.68 for M 2 and S 2 , respectively).Regression-based total dissipation for each constituent (Figure S10 in Supporting Information S1), calculated as the sum of D deep and D shelf from the regression models for each constituent, compares well with explicitly modeled dissipation (r > 0.90) for ICE-6G.A similarly good t can be achieved with the uniform SL simulations.However, it is notable that for M 2 , the uniform SL simulations show only very Interestingly, pronounced differences emerge in the relationship between GMSL and M 2 open ocean dissipation (Figure 8) across the sets of simulations.Whilst both ICE-6G and the uniform SL simulations show a strongly negative relationship, the slope of the relationship is very different, with the uniform SL simulations showing much smaller increases in dissipation than the ICE-6G simulations for the equivalent GMSL decrease.The Gowan simulations fall somewhat between the two other estimates.These differences may reflect the strong control exerted by coastline positions on the large M 2 open ocean tides (e.g., Arbic et al., 2009;Green, 2010;Wilmes et al., 2019), especially as differences in ice sheet extent and non-uniform sea level changes cause offsets in coastlines in comparison to the uniform sea level drop case.On the other hand, the non-uniform sea level changes driven by GIA processes, which are especially pronounced close to ice sheets where we generally see the largest tides, may also be contributing to these differences (see, e.g., Arbic et al., 2008).

Dissipation Estimates Over the Last 430 kyr
Following on, we now use the SL16 sea level curve and associated errors (Figure 9a) to infer dissipation back to 430 ka by applying all three regression models.The SL16 sea level curve shows large variability in GMSL with sea level fluctuating between highstands of ∼ −8 m and +16 m during the interglacials and sea level lowstands of 98-129 m during glacial maxima.On average, GMSL was 54 m lower than present over the last 431 ka, according to SL16.The SL16 GMSL estimate compares well with the ICE-6G derived GMSL for the last glacial cycle (r = 0.98, p = 0.00) (Figure 9a).Small offsets (∼15 m) can be seen during MIS 5, when climate entered into the last glacial phase.However, large differences emerge between SL16 and Gowan, which are especially pronounced during MIS 3 and late MIS 5.
Regression-inferred dissipation values for the last glacial cycle closely follow the explicitly modeled values both for ICE-6G and Gowan (Figures 9b-9k).The largest offsets occur when the sea level curves of ICE-6G and Gowan disagree with SL16, for example, around the relative sea level highstands at 100 ka and 80 ka for ICE-6G, and during MIS 3 for Gowan.The disparities are most pronounced for M 2 open ocean dissipation.Looking back over the last four glacial cycles, the ICE-6G regression model suggests that M 2 open ocean dissipation was strongly increased with respect to present-on average by 120%, that is, more than doubled-apart from during interglacial phases.When GMSL drops by more than 100 m, M 2 open dissipation increases on average by a factor 2.6.For the Gowan model, the mean M 2 open ocean enhancements (36%) are less pronounced and open ocean dissipation during sea level lowstands increased by a factor 1.7.These Gowan change estimates are very similar to those obtained using the uniform SL model (Figure S11 in Supporting Information S1).For shelf dissipation, all three models give very similar results, estimating the mean decrease in dissipation to be around 20% over the last 430 kyr.Total M 2 dissipation averaged over the last 430 kyr was between 3% (uniform SL), 5% (Gowan), and 37% (ICE-6G) greater than at present.For the other constituents, average total dissipation changes are less pronounced, generally within ±25%.However, considerable fluctuations in the partitioning of tidal energy between the open ocean and the shelf seas can be seen; with S 2 , K 1 and O 1 open ocean dissipation increasing during glacial periods and decreasing toward present-day values during interglacials, respectively, and shelf dissipation behaving in an opposing manner.Taken together, these results suggest that globally integrated open ocean dissipation was on average between 28% (Gowan) and 85% (ICE-6G) greater than at present, with enhancements by a factor 1.5 to 2.2 during sea level lowstands and values similar to present during highstands.However, it is worth noting that the exact magnitude of the estimated changes is dependent on the specific sea level (and ice sheet) model adopted.

Discussion and Implications
Our simulations are able to reproduce the findings of other studies for the LGM, deglacial and Holocene (Arbic et al., 2004;Green, 2010;Griths & Peltier, 2008, 2009;Uehara et al., 2006;Wilmes et al., 2019Wilmes et al., , 2022;;Wilmes & Green, 2014), but also provide the first continuous estimates of tides and tidal dissipation for the last 430 kyr, thus expanding our knowledge on tidal dynamics prior to the LGM.
Our results show that, apart from interglacial phases, tides and tidal dissipation were different from present during most of this extended time period, covering multiple glacial-interglacial cycles.Changes were especially pronounced for the M 2 tidal constituent, which displays near-resonant behavior during the LGM (Arbic et al., 2004;Green, 2010;Griths & Peltier, 2008, 2009;Uehara et al., 2006;Wilmes et al., 2019;Wilmes & Green, 2014).Our results show that the M 2 tide in the open ocean was strongly enhanced with respect to present for most of the last glacial cycles, whereas for the constituents S 2 , K 1 , and O 1 the relative enhancements were smaller and conned to periods with the lowest sea levels.Peak open ocean amplitudes and dissipation occurred when GMSL dropped below 70 m.Notably, for the LGM and for both the ICE-6G and Gowan simulations, peak dissipation values do not coincide with the lowest sea levels and greatest ice sheet extent but reach their maximum 4-5 kyr after the lowstand as the ice begins to recede and sea level begins to increase (at 22 ka for ICE-6G and 15 ka for Gowan).This is likely related to the location of grounded ice margins and thus regional shelf area which in turn affects resonance properties of the glacial Atlantic (see experiments in Wilmes et al., 2019, which show that reduced Weddell Sea ice extent increases Atlantic dissipation by ∼1 TW).
However, there are considerable differences between the two reconstructions for the last glacial cycle, and thus differences in the simulated tides.Whilst the two reconstructions show relatively similar GMSL and ice sheet extent from 20 ka to present and agree during the two peak glaciations (LGM and MIS 4), they strongly differ in ice sheet extent and GMSL during MIS 3 (55 ka to 25 ka) and late MIS 5 (80 ka to 65 ka) (see Figure 2).The ICE-6G reconstruction shows an extensive Laurentide ice sheet (LIS) from 80 ka and Fennoscandia is largely covered by ice from 70 ka onwards.In contrast, the Gowan reconstruction has extensive ice sheets over North America and Fennoscandia only during the glacial periods around 60 ka and the LGM.Ice sheet extent leading up to these periods (80-65 ka and 55-30 ka) is much reduced in comparison to ICE-6G, and restricted to smaller more localized ice caps (see Figure 2).This leads to pronounced differences in local relative sea level of over 100 m and GMSL offsets of over 65 m between the two reconstructions during MIS 3 (see Figures 1a and 2).GMSL during MIS three is subject to high uncertainty with sea level estimates for this periods remaining debated (e.g., Dalton et al., 2016Dalton et al., , 2022;;Pico et al., 2016;Pico et al., 2017).Because of a lack of direct markers of GMSL prior to the LGM, pre-LGM GMSL is often inferred from proxy records such as δ 18 O from benthic foraminifera (e.g., Waelbroeck et al., 2002) or planktonic foraminifera (Shakun et al., 2015), from spleotherms in the Red Sea (Grant et al., 2014;Rohling et al., 2009) and compilations from multiple statistically analyzed records (Spratt & Lisiecki, 2016).However, because δ 18 O carries the imprint of both global ice-volume and ocean-temperature changes, the records need to be corrected for temperature changes to obtain sea level curves (e.g., De Boer et al., 2014).Large offsets in MIS 3 GMSL of 30-60 m have been found between ice sheet based reconstructions (e.g., Dalton et al., 2016;Dalton et al., 2022;Gowan et al., 2021;Pico et al., 2016;Pico et al., 2017) and records inferred from marine δ 18 O (e.g., Spratt & Lisiecki, 2016;Waelbroeck et al., 2002), with ice sheet reconstructions suggesting that (a) the Laurentide Ice Sheet during MIS 3 was much less extensive than previously thought and (b) GMSL was ∼30-50 m lower than at present (Dalton et al., 2016(Dalton et al., , 2022;;Pico et al., 2016Pico et al., , 2017)).This implies that, for MIS 3, the Gowan sea level and ice sheet reconstruction may be more appropriate as it captures the relative sea level highstand and reduced ice sheet extent.
By comparing the ICE-6G and Gowan tide simulations for the LGM, it becomes apparent how sensitive the tidal dynamics are to relative small changes in bathymetry, that is, changes in both sea level and ocean basin shape (e.g., through ice sheet extent changes) when in a near-resonant state.At 20 ka, M 2 open ocean dissipation differs by 0.5 TW between the ICE-6G and Gowan simulation, despite a very similar GMSL drop of 116 and 120 m, respectively.Experiments in Wilmes et al. (2019) showed that LGM M 2 tidal dynamics were remarkably insensitive to offsets in GMSL (±10 m) but that changes in the location of the land-ocean boundaries could have dramatic effects (e.g., altering ice sheet extent in the Weddell Sea led to a >1 TW change in M 2 dissipation).This is also shown when comparing the uniform SL change simulations with those that have realistic glacial land ocean boundaries (1.8 TW vs. 2.5 and 2.0 TW for ICE-6G and Gowan, respectively).Comparing the ICE-6G and Gowan 20 ka bathymetries (see Figure 2) shows that the Gowan reconstruction has more ice around the margins of the Laurentide and the Greenland Ice Sheets but also around the margins of Antarctica.Furthermore, whilst GMSL is very similar between the two reconstructions, regional sea levels in the Labrador Sea differ by over 50 m (Figure 2) which may also influence tidal dynamics and contribute to the dissipation offsets.It is also worth noting that the LGM sea level lowstand does not occur at the same time in the two reconstructions: in ICE-6G, the lowest sea levels occur around 26 ka, whereas for Gowan they occur at 20 ka, and the ICE-6G lowstand is 10 m lower than in Gowan (−130 m vs. −120 m).
This work and that of Wilmes et al. (2019) demonstrates that small changes in land-ocean boundaries, for example, through changes in ice sheet extent, could dramatically alter open ocean tidal energy dissipation on the order of ±1 TW.In experiments using the intermediate climate model UVic (Schmittner et al., 2015;Wilmes et al., 2019Wilmes et al., , 2021)), changes of this order of magnitude correspond to a ∼2-4 Sv change in AMOC strength.These strong sensitivities in glacial tides could cause rapid changes in the amount of tidal energy available for ocean mixing which, in turn, could impact ocean circulation and thus climate.Feedbacks between tides, ice sheet extent, ocean mixing, and climate may be a modulating factor for glacial climate, especially for period such as Heinrich events, where ice loss from Hudson Strait has been postulated to have affected tidal dynamics and vice versa (e.g., Arbic et al., 2004;Velay-Vitow et al., 2020).For instance, in the Gowan simulations, between 60 ka and 25 ka, fluctuations in M 2 open ocean dissipation can be seen which appear related to Hudson Bay & Strait dissipation, thus suggesting that Hudson Bay & Strait ice cover may aect (North) Atlantic tidal energy availability.Furthermore, the large tides along the ice sheet margins-both for the Antarctic and Northern Hemisphere ice sheets-may have contributed to the rapid ice sheet retreat linked with the deglacial increase in sea level (Batchelor et al., 2023;Gomez et al., 2020), for example, through ice shelf tidal flexure and crevasse formation (e.g., Olinger et al., 2019).
Because the sea level simulations used for the last glacial cycle include no ice history prior to 122 ka (ICE-6G) and 80 ka (Gowan), the oldest time slices (i.e., LIG for ICE-6G and late MIS 5 for Gowan) are likely associated with the largest errors, as each time step is affected by previous ice sheet and sea level history.For example, because no ice history was present for the penultimate glacial maximum, the exact sea levels and thus tides during the LIG, are likely subject to larger errors.It is also possible that, due to ice history "memory" from previous glaciations, tidal dynamics moving into a glaciation may behave differently than during the deglacial phase.However, at present, due to a lack of global high resolution ice sheet reconstructions beyond the LIG, it is not possible to evaluate this point.Another source of uncertainty is the relatively low temporal resolution of the ice sheet histories (2.5 kyr for Gowan throughout the last glacial cycle and 2 kyr prior to 32 ka for ICE-6G, ≤1 kyr after 32 ka) and for ICE-6G, the horizontal resolution of the ice sheets (1/2° × 1/2°).It is conceivable that the adopted histories lack small scale ice sheet and bathymetry fluctuations, which, during tidally resonant states, could have led to large short-term regional and supraregional shifts in tidal dynamics.
Simulations of climate during the last glacial cycle or parts thereof generally neglect transient changes in tidal energy input to ocean mixing and hold mixing rates constant (e.g., Menviel et al., 2017).However, this work emphasizes that tidal dissipation strongly varied between glacial and interglacial phases and that large short-term variations may have occurred during transitions into and out of glacial phases.The explicitly simulated glacial cycle tidal dissipation time slices can be used to adjust vertical mixing in transient climate model simulations.These also account for regional responses in dissipation which may differ from the global mean.Specifically, during MIS 3, dissipation enhancements in the Southern Hemisphere are greater than in the Northern Hemisphere, which stands in contrast to full glacial conditions during the LGM.

Conclusion
In this study, we have simulated tidal amplitudes and tidal dissipation for the last glacial cycle encompassing the LIG to present.In addition, we have estimated tidal dissipation for the last four glacial cycles spanning back to 430 ka.Our findings suggest that tides likely differed from present during most glacial phases, with maximum tidal amplitudes and open ocean dissipation occurring during sea level lowstands due to the resonant properties of the glacial ocean with regards to the semi-diurnal tidal forcing.Additionally, we discovered that during glacial phases, semi-diurnal tides, particularly M 2 , were highly sensitive to sea level and ice sheet extent, which remain uncertain prior to the LGM.These results are of wider relevance as they indicate that the availability of tidal energy during the late Pleistocene strongly differed from the present, potentially impacting ocean mixing, and thus ocean circulation and climate in the past.
The tidal simulations and relevant grid files together with the Matlab routines necessary to read the data are available at the repository Zenodo (Wilmes et al., 2023a).The dissipation time slices are also accessible via Zenodo (Wilmes et al., 2023b).

Figure 1 .
Figure 1.(a) Global mean sea level (black line) and shelf area (in % of global ocean area) for ICE-6G (solid lines) and Gowan et al. (2021) (dotted lines).Marine isotope stages 1-5 are indicated by shading and labels above the timeline.For MIS 5, letters indicate substages.Interglacial periods are highlighted by darker yellow shading.Italic labels show climate periods: LIG = Last Interglacial; LGM = Last Glacial Maximum.(b-e) Dissipation for the constituents M 2 , S 2 , K 1 and O 1 , respectively.Open ocean (deep) dissipation is plotted in blue, shelf sea dissipation in red and globally integrated dissipation in black.ICE-6 G values are plotted with solid lines and Gowan et al. (2021) values with dashed lines.(f) same as (b-e) but for the sum of the constituents.The dashed thin straight lines give present-day values as a reference.
presentation of the two ice sheet and sea level reconstructions, seeGowan et al. (2021), andArgus et al. (2014),Peltier et al. (2015) andPedersen et al. (in prep).Whilst a fully developed Laurentide Ice Sheet is present from around 100 ka in ICE-6G, a similarly extensive ice sheet is not formed until after 65 ka in the Gowan reconstruction.This leads to a GMSL offset of ∼30 m between the two reconstructions during late MIS 5 and early MIS 4 and large differences in shallow shelf sea area.During the MIS four peak glaciation, the ICE-6G Laurentide

Figure 2 .
Figure 2. Differences in sea level and ice sheet extent between ICE-6G bathymetries and Gowan et al. (2021) for selected time slices.Differences in sea level change plotted as ICE-6G minus Gowan (i.e., blue colors indicate lower sea levels in ICE-6G and red colors higher sea levels in ICE-6G than in Gowan).ICE-6G ice sheet extent is indicated with orange hashing, Gowan ice sheet extent is given with blue hashing.

Figure 4 .
Figure 4. Same as Figure 3 but for the Gowan simulations.

Figure 5 .
Figure 5. (top left panel) Present-day tidal dissipation plotted as the sum of the dissipation of the constituents M 2 , S 2 , K 1 and O 1 .(other panels) Dissipation (as sum of constituents M 2 , S 2 , K 1 , and O 1 ) difference with respect to present for selected time slices over the last glacial cycle for the ICE-6G simulations.Land is shaded in gray.

Figure 6 .
Figure 6.(top left) Present-day tidal dissipation plotted as the sum of the dissipation of the constituentsM 2 , S 2 , K 1 , and O 1 .(otherpanels) Dissipation (as sum of constituents M 2 , S 2 , K 1 , and O 1 ) difference with respect to present for selected time slices over the last glacial cycle for the Gowan simulations.Land is shaded in gray.

Figure 7 .
Figure 7. Correlation coefficients between GMSL and dissipation by tidal constituent from runs with ICE-6G bathymetry (circular markers), with uniform sea level changes (squares) and the Gowan bathymetries (crosses).Blue markers are used for open ocean dissipation, red ones for shelf sea dissipation and black ones for globally integrated dissipation.

Figure 8 .
Figure 8. Open ocean dissipation (blue, left column) and shelf sea dissipation (red, right column) plotted against GMSL for the constituents M 2 , S 2 , K 1 and O 1 (top to bottom).Circular makers show dissipation from the ICE-6G simulations, Gowan simulations are shown with crosses, and squares plot dissipation from the uniform SL runs (M 2 and K 1 only).Regression lines are plotted in black (solid for ICE-6G and dashed for uniform SL runs).The correlation coefficients r between the explicitly modeled and regression-estimated dissipation values are printed in each panel for ICE-6G (r I ), the uniform SL simulations (r U ), and the Gowan simulations (r G ). Stars indicate that the correlation is significant at the 99% confidence level.

Figure 9 .
Figure 9. (a) Global mean sea level from the Spratt and Lisiecki (2016) reconstruction (blue line) with shading indicating upper and lower uncertainties, from the ICE-6G bathymetry (black solid line) and from Gowan et al. (2021) (black dotted line).Light blue and light yellow shading with corresponding numbering indicates Marine Isotope Stages.(b, d, f, h, j) Regression estimated dissipation for the constituents M 2 , S 2 , K 1 and O 1 , respectively, and the sum thereof.Open ocean (deep) dissipation is plotted in blue, shelf sea dissipation in red and globally integrated dissipation in black using the ICE-6G regression models.Shading indicates one standard deviation uncertainty of the dissipation estimates.Explicitly modeled dissipation values for the last glacial cycle are overlain in with dashed lines.Light blue and light yellow shading with corresponding numbering indicates Marine Isotope Stages.(c, e, g, i, j, k) same as (b, d, f, h, j) but for the Gowan et al. (2021) regression model.
Note.Amplitude root-mean square errors (RMSE) against TPXO9 for the deep ocean (h > 500 m), shelf seas (h < 500 m), and the global ocean.Integrated dissipation values (diss.)are also given for the deep ocean (h > 500 m), shelf seas (h < 500 m), and the global ocean for both the present-day control and for TPXO9.