Enhanced vertical mixing in the glacial ocean inferred from sedimentary carbon isotopes

Reconstructing the circulation, mixing and carbon content of the Last Glacial Maximum ocean remains challenging. Recent hypotheses suggest that a shoaled Atlantic meridional overturning circulation or increased stratification would have reduced vertical mixing, isolated the abyssal ocean and increased carbon storage, thus contributing to lower atmospheric CO2 concentrations. Here, using an ensemble of ocean simulations, we evaluate impacts of changes in tidal energy dissipation due to lower sea levels on ocean mixing, circulation, and carbon isotope distributions. We find that increased tidal mixing strengthens deep ocean flow rates and decreases vertical gradients of radiocarbon and δ13C in the deep Atlantic. Simulations with a shallower overturning circulation and more vigorous mixing fit sediment isotope data best. Our results, which are conservative, provide observational support that vertical mixing in the glacial Atlantic may have been enhanced due to more vigorous tidal dissipation, despite shoaling of the overturning circulation and increases in stratification. Vertical mixing in the Atlantic Ocean during the Last Glacial Maximum may have been enhanced, not reduced, due to more vigorous tidal dissipation, suggest an ensemble of ocean simulations and glacial sediment carbon isotope records.

C irculation and mixing of the deep ocean have been suggested as important drivers of glacial-interglacial changes in ocean carbon storage and atmospheric CO 2 . A number of studies hypothesize decreases in vertical mixing during the Last Glacial Maximum  LGM) due to an increase in stratification [1][2][3][4] , or a shoaling of the Atlantic meridional overturning circulation (AMOC) [5][6][7] , which would have isolated the deep ocean, increased carbon storage there and thus contributed to lowering atmospheric CO 2 . Other studies, however, have suggested enhanced tidal energy dissipation in the deep ocean, which would have increased mixing [8][9][10][11][12][13][14][15] , but the effects of changes in tidal mixing on ocean biochemical cycles and isotope distributions remain unknown and data constraints on deep LGM ocean mixing unexplored.
Tides provide an important energy source for internal mixing that sustains the ocean's meridional overturning circulation 16,17 . Surface (barotropic) tides lose energy through interactions with topographic features on the sea floor. A portion of that energy is dissipated locally, leading to turbulence and mixing there (nearfield tidal mixing) 18,19 . The remainder of the energy is converted to internal waves (far-field internal tides) that propagate away from the generation sites and lead to mixing and dissipation elsewhere 19 . Since tidal mixing arises through interactions with topography, it is greatest near the sea floor and decreases towards the surface. It has been suggested 6 that a shallower interface between North Atlantic Deep Water (NADW) and Antarctic Bottom Water (AABW), as observed for the LGM, reduced mixing between the two water masses, thus isolating the deep ocean and increasing its carbon content. However, LGM sea level was approximately 120-130 m lower than at present and large additional ice sheets covered North America and northern Europe. The resulting changes in ocean depth and basin shape enhanced open-ocean tidal amplitudes, especially for the dominant semidiurnal tidal constituent M 2 , and increased global tidal energy dissipation in the deep ocean during the LGM, with the most prominent enhancements in the North and South Atlantic. At present, the magnitude of the increase is not well constrained and highly sensitive to uncertainties in reconstructions of glacial ice sheet extent 13,14 . Stronger tidal mixing and increased diapycnal diffusion could have strengthened the meridional overturning circulation 14,20 .
While there is agreement that the LGM AMOC was shallower than today [21][22][23][24] , its strength remains a contentious subject with estimates ranging from a weaker 7,22,24,25 to a similar or stronger circulation 23,26,27 . Most studies assessing model versus paleo data to reconstruct the meridional overturning circulation have not considered tidal mixing changes 25,28,29 , or if they did, simplified horizontally nonvarying diffusivities were used 23,30 .
The aims of this study are to (1) evaluate the possible impact of LGM near-field tidal mixing changes on carbon isotope distributions in the glacial Atlantic, (2) to assess whether the existing carbon isotope measurements from LGM sediments provide constraints on vertical mixing, and (3) consider uncertainties in AMOC and tidal energy dissipation. Building on previous work 14, 24 and using an intermediate-complexity global coupled physics and biogeochemistry model that includes near-field tidal mixing, we create an ensemble of LGM circulation states forced with different tidal dissipation estimates and modified atmospheric meridional moisture transport. The latter alters the buoyancy forcing in the Southern Ocean. Since changes in farfield tidal dissipation are not included in our model, the results are conservative. The simulations are compared against radiocarbon and δ 13 C distributions from LGM sediments 7,24,[31][32][33][34] . While these data constrain AMOC depth efficiently they do not provide good constraints on AMOC strength 35 . Our results show that simulations with a shallower and weaker AMOC and more vigorous tidal mixing fit the LGM carbon isotope data best. This suggests that vertical mixing in the glacial Atlantic Ocean may have been enhanced, and not reduced.

Results
Near-field tidal mixing effects on Atlantic circulation and water mass structure. Sampling AMOC and tide uncertainties, we have generated an ensemble of LGM model circulation states by modifying the Southern Hemisphere moisture diffusivity (μ SH ) 24 and applying three different tidal dissipation cases from Wilmes et al. 14 : present day (PD) with 1 TW (1 TW = 10 12 W) global internal tide dissipation, LGM ICE-6G (1.8 TW) and LGM ICE-5G (3.6 TW) (see "Methods" for details and Fig. 1 for run names). The two LGM tide cases result mainly from differences in land ice extent reconstructions, such as in the Weddell Sea, where the position of the grounding line remains uncertain. Figure 1 shows Atlantic streamfunctions for the ensemble. For PD tide forcing, AMOC strengths (here, maximum northward streamfunction at 25°N) range from a shutdown circulation state to an AMOC strength of 13.2 Sv (1 Sv = 10 6 m 3 s −1 ). For the simulations with LGM tide forcing, AMOC strengths range from 7.4 Sv to 15.1 Sv (ICE-6G tides) and 8.8 Sv to 16.7 Sv (ICE-5G tides). Including LGM tidal mixing increases AMOC strengths by around 1-3 Sv. For equivalent AMOC strengths, the AABW circulation cell is enhanced in simulations with stronger tidal forcing (by~0.5 Sv for ICE-6G tides and~1 Sv for ICE-5G tides; see Fig. 2a). This leads to a shoaling of the AMOC by between 200 and 300 m in the cases with stronger tidal forcing in comparison with an equivalent-strength AMOC model forced with PD tides (Fig. 2b). Nevertheless, AMOC strength and depth are highly correlated in our LGM experiments, which makes it difficult to distinguish between the two properties.
Mid-depth (2000-4000 m) increases in stratification, expressed as the buoyancy frequency squared (N 2 ), exceeding 50% relative to the preindustrial control (PIC) are seen for the simulations with AMOC strengths below 11 Sv and corresponding AMOC depths of less than 3000 m in all ocean basins (Figs. 3c and d and Fig. 2). Stratification also affects modeled diffusivities (see equations (2) and (1), respectively). This can be appreciated by comparing simulations with different AMOC states but identical tidal energy dissipation ( Fig. 3c and d). Simulations with a strong and deep AMOC (e.g., i5gT 17Sv) show reduced stratification (N 2 ) and higher k v compared with weak and shallow AMOC cases (e.g., i5g 10Sv). However, the effects of stratification on k v are generally smaller than the effects of tidal energy dissipation ( Fig. 3e and f), such that all simulations with LGM tides have higher diffusivities in the deep Atlantic than any of the simulations with PD tides. In the simulations using LGM tidal mixing, the relative increases in the energy-dissipation rate (ϵ) are much greater (up to 100% and 400% for ICE-6G globally and for the Atlantic, respectively, and more than twice that for ICE-5G tides), thus outcompeting the changes in N 2 and increasing k v at all depths (Fig. 3).
For most simulations with fixed tidal energy dissipation, diffusivities at the interface between NADW and AABW (AMOC depth) decrease as the interface shoals (Fig. 2c), consistent with the hypothesis of Ferrari et al. 6 . However, interface diffusivities are larger for runs with LGM tides than for PD tides at similar AMOC strengths, which can compensate for the effects of AMOC shoaling. This means that, e.g., in the simulations with AMOC strengths of 9 Sv, k v at the interface is enhanced by a factor of 2.5 for ICE-6G tides and 3.4 for ICE-5G tidal mixing.
Model-sediment isotope comparison. Atlantic radiocarbon distributions display strong circulation imprints ( Supplementary  Fig. S4). Generally, northward of~30 ∘ S, southern-sourced waters occupying the lower part of the water column have older radiocarbon ages than waters sourced from the north. Southward of 30 ∘ S, the entire water column is dominated by southern-sourced waters, and its structure, in contrast to the north, shows little change in response to the circulation changes ( Supplementary  Fig. S8). Weak and shallow AMOC cases have older radiocarbon ages in the deep ocean than those with stronger and deeper AMOC. Simulations with LGM tidal mixing show more gradual gradients separating southern-sourced from northern-sourced water masses in the central and North Atlantic compared with PD tide models, for which gradients are sharper. Water masses occupying the deep North Atlantic have younger ages in the LGM tidal mixing cases than equivalent AMOC strength models with PD tidal mixing. A very similar picture emerges for the δ 13 C distributions (Supplementary Figs. S5 and S7). Again, a strong circulation imprint is expressed, and cases with more mixing show reduced vertical gradients and less negative δ 13 C values at depth for an equivalent AMOC strength. The largest discrepancies between modeled and sediment isotope compositions are seen in the deep South Atlantic where modeled values are isotopically too light and ages too young.
Calculations of root-mean-square (RMS) errors between modeled and reconstructed radiocarbon ages and δ 13 C reveal a strong dependence on AMOC strength/depth, both globally and for the Atlantic (Fig. 4a-d). The lowest RMS errors are found for AMOC strengths of 8-10 Sv and AMOC depths of 2.4-2.7 km, and the highest errors for the cases with either a shutdown or a strong and deep AMOC. Simulations with stronger tidal mixing consistently have lower δ 13 C RMS errors in the best-fitting AMOC strength/depth range than the simulations with PD tides (Fig. 4b, c and d). Simulations i5gT 9Sv, i5gT 10Sv and i6gT 9Sv show the best fit in the Atlantic for both radiocarbon ages and δ 13 C. Simulation i5gT 9Sv decreases the global δ 13 C RMS error by 8%, and in the Atlantic between 30°S and 50°N it improves the fit with the radiocarbon and δ 13 C data by 15% and 9%, respectively, in comparison with the best-fitting run with PD tidal mixing and an AMOC strength of 9 Sv. We argue that Atlantic data provide more direct constraints on these simulations than global data since the largest changes in mixing occur in the Atlantic, and Pacific isotope distributions will be affected by additional uncertainties such as changes in deep/intermediate water formation in the North Pacific 36 that were not sampled in our ensemble. More data are needed to reliably represent the deep South Atlantic and Atlantic Southern Ocean where data coverage is sparse and reconstructions of δ 13 C are hampered by higher errors for the present 37 . It is further possible that some climate processes in the South Atlantic and Atlantic Southern Ocean are not fully represented as the very old radiocarbon and negative δ 13 C values in the LGM observations are not reached in this ensemble.
The depth profiles in Fig. 4e-f and Supplementary Figs. S6-S9 illustrate a transition zone between~2000 and 3500 m northward of the equator in which the water mass characteristics change from northern sourced to southern sourced. The vertical gradient in this zone is quite smooth in the observations. It becomes evident that models with weak AMOC states have too strong mid-depth gradients (e.g., pdT 7Sv or i6G 7Sv), whereas simulations with an AMOC strength >12 Sv show a gradient that is too weak, and the profiles have too young radiocarbon ages and too positive δ 13 C values at depth. Only three simulations, those with AMOC strengths of 8.5-10 Sv and ICE-5G and ICE-6G tidal forcing, lie fully within a standard deviation of the mean radiocarbon and δ 13 C depth profiles (see Supplementary Fig. S6). They feature AMOC depths of~2400-2600 m and Atlantic vertical diffusivities that are enhanced by a factor of 2-3.1 in comparison with the simulation with equivalent AMOC depth but present-day tidal forcing. The simulation with ICE-5G tidal mixing and an AMOC strength of 9 Sv shows the best fit for both carbon isotopes, giving a consistently close fit to the mean depth profiles together with the lowest overall RMS errors, both in the Atlantic and globally.
The fit between the LGM simulations and the carbon isotopes is largely determined by AMOC strength/depth, as indicated by the near collapse of all data onto roughly one curve in Fig. 4a-ddespite the large differences in tidal mixing. However, mixing effects emerge for experiments with best-fitting AMOC strengths/ depths for which LGM tides improve the agreement with the carbon isotopes compared with PD tides. This indicates that tracer distributions are influenced by vertical mixing, even if that influence is more subtle than the dominating effect of AMOC structure and thus advection.

Discussion
Shoaling of the NADW-AABW interface moves it away from the high-mixing regions near the sea floor and thus decreases diffusive exchange between those water masses. This effect 6 is captured in our experiments that show a general trend toward lower diffusivities with shallower AMOC depths for simulations with constant tidal dissipation (Fig. 2c). However, increased LGM tidal energy dissipation over-compensates for this effect, such that the three best-fitting LGM simulations (i5gT 9Sv, i5gT 10Sv, and i6gT 9Sv) have higher diffusivities on the NADW-AABW interface compared with the PIC simulation despite a shallower interface. In these runs, the remotely dissipated component is kept constant, however, it would be expected to increase in line with the locally dissipated tidal energy since most of the energy from the barotropic tide is transferred to internal waves that propagate away from the generation sites 19 . Considering these far-field effects would thus further increase LGM diffusivities. Our model assumes that far-field effects account for two-thirds of the total tidal energy dissipation. Although the partitioning into near-and far fields may be different 18 , it would suggest that the effects of increased tidal mixing could be a factor of two larger than simulated here. On the other hand, the model underestimates present-day AMOC depth and thus the degree of shoaling during the LGM. Assuming a modern AMOC depth at 25 ∘ N of 4.4 km from observations ( Supplementary Fig. S1 38 ), the PIC model predicts higher interface diffusivities in the modern tropical (30°S-30°N) Atlantic by about a factor of 2. Including higher latitudes leads to a smaller difference.
Our model-data comparison implies that vertical mixing rates in the glacial Atlantic, including the Atlantic sector of the Southern Ocean ( Supplementary Fig. S10e), may have been enhanced. This contradicts conclusions from a number of studies that suggest reduced vertical mixing due to increased density stratification was an important factor in lowering atmospheric CO 2 during the LGM 1,6,2-4,39 . While our simulations do not infer direct mechanisms for the glacial atmospheric CO 2 drop, they include the effect of stratification changes on diffusivities and show that those effects could have been dwarfed by the impacts of increased tidal energy dissipation, which were not included by these studies. This therefore implies that reduced Atlantic vertical diffusivities may not have been the mechanism for increased Atlantic carbon storage, rather that carbon storage was enhanced, despite of more vigorous vertical mixing. Moreover, previous work 1 used potential density referenced to the surface to calculate changes in N 2 , which overestimates the effects on k v compared with using locally referenced potential density 40 . Based on analysis of oxygen isotope data, it has been suggested that the ratio between lateral transport and vertical mixing (Ψ/κ) during the LGM was increased by at least a factor of two compared with present day 5 . Here, we show that when we consider tidal mixing consistent with LGM ocean basin shape and depth, global diffusivities increase by a factor of 2-4 throughout the water column, with the strongest enhancements between 2000 and 4000 m. Compared with the modern ocean, the AMOC is reduced by about half in our best fitting simulations, which suggests Ψ/κ was lower by a factor of 4-8, in contrast to Lund et al.'s 5 results. Our model does not include oxygen isotopes and thus we cannot assess the consistency of our simulations with those data.
More work is needed to quantify the effect of tidal mixing changes on carbon storage and atmospheric CO 2 . Our best-fitting runs show similar carbon-storage values as the best fitting state in Muglia et al. 24 (within~120 GtC; Supplementary Fig. S9), which was shown to explain more than three-quarters of the glacial-interglacial CO 2 difference 41 . The same method could be used to analyze the experiments presented here, but based on the similarity in tracer distributions, we would expect similar effects on CO 2 .
The AMOC strength may not be constrained well by the data 35 such that a shallow and strong AMOC state, which was not tested here, is possible 23 . However, compared with a shallow and weak AMOC state, a shallow and strong circulation would sharpen the mid-depth isotope gradients, which would require even more vertical mixing in order to fit the smooth gradient observed in the sediment data. It has been shown 42 that global ocean carbon content is more sensitive to AMOC depth than to AMOC strength. Thus, a strong and shallow AMOC with increased vertical mixing would presumably have a similar carbon content than the weak and shallow cases considered here.
It has been suggested 43,44 , using a highly simplified ocean model setup with uniform globally uniform k v , that in the present day, when diapycnal mixing rates are low, the mid-depth interhemispheric circulation is primarily set by the winds over the Southern Ocean and the associated upwelling they drive, and that the circulation structure is relatively insensitive to small changes in diapycnal mixing. However, it has been demonstrated 43,44 that a second circulation state may exist in the presence of much higher diapycnal mixing rates where the circulation and stratification structure are set by the diapycnal diffusivity and North Atlantic surface buoyancy distribution. The effect of the increased diapycnal mixing on the isopycnals (Fig. 9 in Nikurashin and Vallis 44 ) can be seen when simulations pdT 9Sv and i5gT 9Sv are compared (Fig. S11), suggesting that altered LGM tides may be shifting the circulation from primarily wind driven to more diapynal mixing driven during the LGM.
The effects of increased tidal energy dissipation are likely underestimated in our model because we only consider the locally dissipated energy, while remotely dissipated energy is kept constant in the term k bg . There are uncertainties in the proportions of near-field and far-field energy dissipation, and recent research points toward a spatially varying partitioning 18 . Furthermore, our current setup uses a fixed-exponential vertical decay scale, however, dynamic decay scales depending, e.g., on stratification have been proposed 45,46 . Future research should implement a parameterization of far-field effects 47 , dynamic partitioning of locally and remotely dissipated energy, and allow the vertical decay scale to evolve between climates.
For the best-fitting circulation cases, increased LGM tidal mixing improves the fit to the Atlantic radiocarbon and δ 13 C sediment data by around~10-15%. This implies, contrary to the prevailing view in the literature [1][2][3][4][5][6]25 that LGM Atlantic vertical mixing may have been strengthened, rather than weakened, and it further highlights that explicit and relevant tidal energy sources need to be implemented in ocean circulation models. We show, for the first time to our knowledge, that sediment data provide constraints on past ocean mixing. The reconstruction method could be applied to other, relatively data-rich paleo periods.

Methods
Climate model. Here, the intermediate-complexity University of Victoria (UVic) 48,49 climate model, version 2.9, is used. The three-dimensional dynamical ocean model has a horizontal resolution of 3.6°× 1.8°and 19 vertical levels, and is coupled to a single-level 2D energy-moisture balance atmospheric component and a thermodynamic sea-ice model that allows sea ice to evolve in space and time. The climate model is run in the same setup as in Muglia et al. 24 , apart from the tidal mixing, and thus includes LGM wind speed anomalies from simulations with fully coupled ocean-atmosphere general circulation models as part of the Paleoclimate Model Intercomparison Project (PMIP) version 3 ensemble.
UVic is coupled to the Model of Ocean Biogeochemistry and Isotopes (MOBI), version 1.8 24,[50][51][52] , which includes prognostic equations for phosphate, nitrate, dissolved and particulate iron, O 2 , dissolved inorganic carbon, dissolved organic and particulate matter, zooplankton, and diazotrophs. MOBI includes an advection-diffusion scheme, so tracers are influenced by both physical ocean dynamics and biogeochemical processes. For details on the prognostic 14 C cycle in the model, see Muglia et al. 24 . We use the same prescribed iron fluxes as Muglia et al. 24 for the preindustrial control simulation. For the LGM, we use their iron fluxes with enhancements by a factor of 10 in the Southern Ocean as these give the best fits to sediment nitrogen and carbon isotope data 24 .
Near-field tidal mixing is implemented following Schmittner and Egbert 53 and as applied in Schmittner et al. 20 and Wilmes et al. 14 . It accounts for subgrid-scale bathymetric effects on the depth and spatial structure of the tidal energy input, and distinguishes between semidiurnal and diurnal tidal constituents. Three different tidal dissipation fields, taken from Wilmes et al. 14 , were used for the LGM climate model simulations: 1. Present-day tidal dissipation (PD) 2.
LGM ICE-6G tidal dissipation (LGM ICE-6G) PD tidal dissipation was calculated from a tide model simulation using present-day ocean bathymetry. The LGM dissipation fields came from tide model runs using two different LGM ocean bathymetries based on the sea-level changes in the ICE-5G (VM L90) version 1.2 54 and ICE-6G (VM5a) 55,56 models. The tidal dissipation fields 14 have a horizontal resolution of 1/8°× 1/8°and include tidal energy losses for the principal diurnal and semi-diurnal tidal constituents M 2 , S 2 , K 1 , and O 1 . These fields were mapped onto the climate model grid using the methodology described in Schmittner and Egbert 53 . In the climate model, the diapycnal diffusivity k v is given by where k bg is the background diffusivity with a value of 0.3 × 10 −4 m 2 s −1 , which includes remotely dissipated tidal energy (far-field) and mixing through other processes. Γ, the mixing efficiency, is set to 0.2, and N 2 is the buoyancy frequency. N 2 is expressed as where g is the acceleration due to gravity, ρ is density, and z is depth in the water column. For details on the calculation of N 2 see Supplementary Methods S3. The rate of tidal energy dissipation, ϵ is given by where D IT;TC ðx; y; z 0 Þ is the energy flux from the barotropic to the internal tide from the high-resolution tide model, D IT , mapped onto the climate model grid accounting for subgrid-scale bathymetric effects in the vertical (thus the dependence on z 0 ). F is the vertical decay function using an e-folding depth of 500 m above the sea floor H. q TC , the local dissipation efficiency, accounts for the critical latitude y c of diurnal and semidiurnal tidal constituents (TC): q TC ¼ 1; for jyj > y c;TC 0:33; otherwise: ð4Þ y c is 30 ∘ for the diurnal constituents K 1 and O 1 , and 72°for the semidiurnal constituents M 2 and S 2 . It is important to note that the tidal mixing parameterization used here considers only changes in the locally dissipated energy, which amounts to only about one-third of the total energy dissipation. The remaining two-thirds is assumed to propagate away as internal waves and dissipate elsewhere. However, we do not account for changes in the remotely dissipated fraction, which is represented by a constant background diffusivity in our simulations. This is important for the interpretation of our results because it will make our simulations underestimate the true effects and thus be conservative. Following Muglia et al. 24 , we stepwise modify the Southern Hemisphere meridional moisture transport F q,SH , which has been shown to alter Antarctic Bottom Water formation and AMOC strength 57 , thus giving a spread of different AMOC states. In the calculation of specific moisture in UVic, the moisture eddy diffusivity term μ accounts for processes not resolved by advection. μ is the largest in the mid-latitudes where eddies play an important role. An offset term μ SH is added in the Southern Hemisphere in order to increase the fit with observations. In our runs, we multiply μ SH by factors of 0, 0.1, 0.25, 0.5, and 1 in order to achieve a suite of different ocean circulation states. Values lower than one decrease the meridional moisture flux and precipitation at high latitudes, which increases Southern Hemisphere sea-ice extent, increases salinity, surface pCO 2 , and carbon content around Antarctica and thus of AABW, and enhances AABW flow rate and decreases the AMOC (see Supplementary Fig. S10a-d).
Aside from the preindustrial control (PIC), we run 15 LGM simulations (for setup details see Supplementary Methods S1), using each of the three tidal forcings (PD, LGM ICE-5G and LGM ICE-6G) in combination with the 5 μ SH offset factors described above.
The PIC is evaluated in Supplementary Methods S2 and Tables S1 and S2. The PIC reproduces observed modern distributions of both isotopes well, such as lower (older) δ 13 C (radiocarbon ages) in southern-sourced and higher (younger) values in northern-sourced deep water (Supplementary Fig. S2). It is also consistent with observational estimates of circulation indices (Supplementary Table S1), although similar to other models, its modern AMOC depth at 25 ∘ N is too shallow (~3.1 km) compared with observational estimates (~4.4 km; see Supplementary Fig. S1 and Fig. 5 in Danabasoglu et al. 58 ), which leads to too low (old) δ 13 C (radiocarbon ages) in the abyssal North Atlantic (Supplementary Fig. S2).
Proxy data. We compare our modeled radiocarbon ages and δ 13 C distributions to present-day and LGM radiocarbon and sediment isotope ratios. In total, 256 LGM radiocarbon age data points come from sediment data compiled by Skinner et al. 7 (see their work for details on the data). Outlier data points with ages greater than 6000 14 Cyr were excluded, leaving a remainder of 246 points. Preindustrial water column radiocarbon ages were taken from the GLODAP 59 database for the locations of the LGM data. In total 434 preindustrial (core-top) and LGM δ 13 C values from Cibicides foraminifera come from a global database 31 . Nine extra data points were added in the South Pacific 33,34 and one in the equatorial Atlantic 32 . In this work, we do not present results on δ 15 N because several studies 24,51 demonstrate that δ 15 N distributions are dominated by biological nutrient-utilization effects rather than circulation effects. An analysis (not shown) comparing LGM minus PIC δ 15 N values for the simulations presented in this study with the δ 15 N data from Muglia et al. 24 gives similar results to Muglia et al. 24 who show that simulations with the same Southern Ocean iron fluxes have similar errors, regardless of circulation structure.

Code availability
The UVic model code and run-setup files, together with the tidal dissipation and μ SH files, are available in Supplementary Data S1.