Mediterranean heat injection to the North Atlantic delayed the intensification of Northern Hemisphere glaciations

The intensification of Northern Hemisphere glaciations at the end of the Pliocene epoch marks one of the most substantial climatic shifts of the Cenozoic. Despite global cooling, sea surface temperatures in the high latitude North Atlantic Ocean rose between 2.9–2.7 million years ago. Here we present sedimentary geochemical proxy data from the Gulf of Cadiz to reconstruct the variability of Mediterranean Outflow Water, an important heat source to the North Atlantic. We find evidence for enhanced production of Mediterranean Outflow from the mid-Pliocene to the late Pliocene which we infer could have driven a sub-surface heat channel into the high-latitude North Atlantic. We then use Earth System Models to constrain the impact of enhanced Mediterranean Outflow production on the northward heat transport in the North Atlantic. In accord with the proxy data, the numerical model results support the formation of a sub-surface channel that pumped heat from the subtropics into the high latitude North Atlantic. We further suggest that this mechanism could have delayed ice sheet growth at the end of the Pliocene. Subsurface heat transfer from the subtropical Atlantic Ocean to high latitudes during the mid-to-late Pliocene may have been responsible for a delay in the intensification of Northern Hemisphere glaciations, according to sedimentary proxy records and numerical modelling.

T he intensified growth of Northern Hemisphere ice sheets (iNHG; ∼3. 2-2.4 Ma (million years ago)) towards the end of the Pliocene marks the onset of the profound glacial-interglacial fluctuations that are characteristic for the subsequent Pleistocene epoch 1

. The acceleration of Northern
Hemisphere ice sheet growth is proposed to have been triggered by a simultaneous global atmospheric cooling trend linked to, firstly, a steady decline in atmospheric carbon dioxide concentration starting at~4.0 Ma (pCO 2 ; see Fig. 1) [2][3][4] . Secondly, it has been proposed that the strength of the global oceanic conveyor belt weakened starting at~4-3 Ma, which lowered the heat and moisture exchange between ocean and atmosphere [5][6][7] . Despite the proposed decline of the oceanic heat transport and pCO 2 during the Plio-Pleistocene transition, which provides ideal background conditions to support the iNHG, the high-latitude North Atlantic (>50°N) experienced a strong sea surface warming (2-3°C relative to present day) between~2.9 and 2.7 Ma (cf. ODP Site 982 in Fig. 1; for other high-latitude SST records see Fig. S1) 8 . This anomalous sea surface temperature (SST) increase most likely led to higher Arctic and sub-Arctic air temperatures, which in turn could have delayed the formation of major ice sheets in the Northern Hemisphere 5,9 . In fact, the intensification of ice sheet development on Greenland, Scandinavia and North America only accelerated after~2.7 Ma during glacials G6-G4 subsequent to the North Atlantic SST warming event ( Fig. 1) 10,11 .
The anomalous surface warming of the high-latitude North Atlantic contrasts the relatively steady decline in pCO 2 betweeñ 3.2 and 2.4 Ma, and thus argues against pCO 2 as a dominant driver of high-latitude SST variability (Fig. 1). The SST increase at high latitudes in the North Atlantic could also result from a transitory invigoration of the Atlantic Meridional Overturning Circulation (AMOC) as its strength modulates the northward directed heat transport 12 . AMOC-driven SST anomalies are expected to lead to synchronous SST changes in mid-(~40°N) and high-latitude (~60°N) regions of the North Atlantic 13 . However, the mid-latitude North Atlantic did not experience a similar amplitude surface warming as the high latitudes (cf. IODP Site U1313 in Fig. 1; for other mid-latitude SST records see Fig. S1) between~2.9 and 2.7 Ma, arguing against a strong oceanic overturning circulation being responsible for the entire magnitude of the reconstructed high-latitude SST warming. Similarly, a potential northward shift of the arctic front, facilitating protrusion of warm subtropical waters into the high latitudes, is at odds with the inferred southward shift of the North Atlantic frontal systems during the Late Pliocene glacials 14 . Alternatively, the high-latitude SST warming between~2.9 and 2.7 Ma could potentially be related to the closure of the Central American Seaway (CAS) and Arctic gateways (e.g. Bering Strait and Canadian Archipelago), which have been argued to facilitate high-latitude warming in the North Atlantic via an intensified Gulfstream/North Atlantic Current 15,16 . However, the anomalous SST warming between~2.9 and 2.7 Ma does not coincide with a presumed closure of the Arctic gateways during the M2 event (~3.3 Ma) 17 and clearly postdates the main phase of constriction of the CAS that occurred from 4.8 to 3.8 Ma 5 . In fact, during the Late Pliocene, and our studied time frame, the CAS sill depth fell most likely below 100 m 5 prior to its final closure, and thus did no longer have a substantial pull on AMOC strength as it presumably did during the Early Pliocene 18,19 . Hence, the observed Fig. 1 Overview of global climate dynamics during the Plio-Pleistocene transition. a Ice volume reconstruction with Glacials G4/6 and G10 denoted 44 (b) different atmospheric pCO 2 reconstructions for the Plio-Pleistocene transition (white circles; 2 orange circles; 3 blue circles 4 ); with LOESS smoother (black line; smoothing factor 0.07) and 95% confidence interval (green shading) derived from a Monte Carlo Simulation with 1000 iterations. c Sea surface temperature (SST) records of mid-latitude Site U1313 (dark red line; thick line indicates 0.5 LOESS smoother; see Fig. 2 for site location) 43 and high-latitude Site ODP 982 (blue line; thick line indicates 0.3 LOESS smoother; see Fig. 2 for site location) 64 on the revised age model of 41 . The grey shading indicates the interval of anomalous SST warming in the high-latitude North Atlantic as depicted by blue arrows. delay in Northern Hemisphere ice sheet growth against the backdrop of global cooling and reduced northward heat transport requires an additional but yet unknown process that specifically alters surface-water temperatures in the high-latitude North Atlantic.
A possible candidate to explain the SST warming of the highlatitude North Atlantic between~2.9 and 2.7 Ma could be the intrusion of dense and warm Mediterranean Outflow Water (MOW) into the central North Atlantic. MOW originates from warm and saline intermediate-to deep-water masses of the eastern and western Mediterranean Sea and enters the North Atlantic via the Strait of Gibraltar 20 at water depths of 800-1200 m (Fig. 2) 21 . Modern observational data show that the injection of Mediterranean-sourced warm and saline water masses into the subtropical intermediate North Atlantic can modify its deep thermocline at~500-800 m (Fig. 2) 22 . Increased MOW presence thereby weakens the stratification at mid-depth in the North Atlantic (Fig. 2) and increases the heat advection along isopycnals from the mid-latitude subsurface to the surface of the subpolar North Atlantic 23 . The formation of this subsurface channel allows heat accumulated at mid-depth of the equatorial and subtropical Atlantic to be transported northwards independent of temperature changes at surface-water level 23 . Hence, increased MOW production and its potential to channel heat at subsurface levels towards the North could explain the anomalous surface warming in the high-latitude North Atlantic during the Plio-Pleistocene transition. Notably, this potential synergetic link between MOW and high-latitude SST variability in the North Atlantic has thus far not been considered in numerical models because they lack the proper representation of the outflow due to the narrow width and shallow depth of the Strait of Gibraltar.
Here we investigate the relationship between MOW production and high-latitude North Atlantic SST variability during the iNHG. We integrate numerical modelling with proxy data to approximate MOW strength and its impact on temperature in the North Atlantic and compare our results with existing SST records. The MOW strength reconstruction is based on X-ray fluorescence (XRF) scanning results from Site U1389 24 located in the Gulf of Cadiz and existing records of hydroclimate variability in the Mediterranean realm. The impact of MOW intrusion on the North Atlantic hydrography is estimated via the atmosphere-ocean-sea-ice-vegetation model COSMOS (Community Earth System Models) under Pliocene boundary conditions.
Results and discussion MOW variability during the Plio-Pleistocene transition. To reconstruct MOW strength, we approximate its flow speed that is proportional to its volume flux 25 by using the sedimentary ln(Zr/ Al) ratio determined by XRF core scanning (see "Methods") 26,27 . This proxy follows the rationale that high bottom flow speeds favour the accumulation of heavy minerals relative to the less dense alumino-silicates 26 . Hence, a high ln(Zr/Al) ratio reflecting a high MOW flow speed also indicates an increased volume flux as well as strength, and vice versa. The temporal ln(Zr/Al) evolution at Site U1389 clearly depicts cyclical changes of the MOW production between 3.1 and 2.55 Ma (see Fig. 3). In fact, the spectral analyses of the ln(Zr/Al) record reveal a persistent precession pacing with a spectral peak at~20 kyr ( Fig. S2; Supplementary Methods). The distinct modulations of MOW strength by precession derive from the well-documented intimate link between MOW and monsoon controlled freshwater discharge into its eastern Mediterranean source region [28][29][30] . There, surfacewater freshening disrupts intermediate convection-and thus the formation of MOW source waters-when monsoonal runoff enters the eastern Mediterranean Sea via the River Nile 28,30 . This situation occurs during precession minima. In contrast, during precession maxima (i.e. low NH summer insolation) the North East African monsoonal system is positioned further to the south, which reduces Nile River runoff and by extension fosters MOW formation 31,32 .
The interval between~2.9 and 2.7 Ma ( Fig. 3b; glacials G14-G7) of Site U1389 is characterized by the repeated occurrence of sandy layers (see Fig. S3). The high sand contents of the drift deposits are particular susceptible of being washed out during the rotary drilling process 33 and prohibited continuous XRF-scanning (gaps in the ln(Zr/Al) record). The coarsening of the sediment between~2.9 and 2.7 Ma implies an intensified MOW production that caused an exceptionally high flow speed regime along the continental slope of the Gulf of Cadiz. This led to increased winnowing and enhanced sand-sized particle deposition as observed at Site U1389 25,26,28,30 . Widespread erosive features and enhanced sand contents testifying to enhanced MOW flow speed during this time interval are also found in seismic profiles and sediment cores across the Gulf of Cadiz 33,34 . Thus, the improved core recovery of the time intervals preceding, and subsequent to, the intensified MOW interval (here termed "iMOW" interval), as evident from the recorded ln(Zr/Al) variability, points to overall reduced sand contents in these drift deposits, and in extension to reduced overall MOW production.
To support our interference of increased MOW production during the iMOW interval, we assessed the fresh water input to its Mediterranean source region by comparing our MOW flow reconstruction to the North African aridity index of ODP Site 967 ( Fig. 3a, b) located in the Levantine Basin ( Fig. 2) 35,36 . We find that the proposed tight coupling between MOW production and North African hydroclimate is supported by the significant correlation between the ln(Zr/Al) record and the North African aridity index for our study interval ( Fig. S4; Supplementary Methods). Interestingly, during the iMOW interval between~2.9 and 2.7 Ma North Africa experienced a strong increase in aridity ( Fig. 3a, b), which is potentially linked to the simultaneous~2.4 Ma eccentricity minimum. In fact, results from numeric modelling as well as proxy data analyses have demonstrated a tight but non-linear linkage between the Mediterranean freshwater budget and changes in Earth's orbital configuration 31,[37][38][39] . During a pronounced eccentricity minimum, the precession amplitude is strongly muted which, according to numerical model and proxy data, causes a substantial weakening of the African monsoonal system 31,40 . This can lead to a reduction of continental runoff into the Mediterranean of up to −2 × 10 15 lyr −1 and thus facilitate an increase of Mediterranean salinity by up to~10 PSU 38 . In fact, we find a maximum~2.5‰ increase of benthic δ 18 O ivf-sw , reflecting intermediate water depth salinities in the western Mediterranean basin, coinciding with the iMOW phase (Fig. 3d) 41 . This indicates that the salinity of the MOW was enhanced by up to~7 PSU during the iMOW phase relative to the time intervals 3.1-2.9 and 2.7-2.6 Ma, respectively (see "Methods" for discussion about the conversion of δ 18 O ivf-sw into salinity) 41 . This supports our finding that during the iMOW phase strong North African aridity facilitated a substantial increase in Mediterranean water salinity and consequently enhanced MOW production as documented by the erosional features along the continual slope of the Gulf of Cadiz.
Influence of MOW on the North Atlantic heat distribution. To assess the potential effect of increased MOW intrusion on the spatial heat distribution of the mid-latitude North Atlantic, and its potential for explaining the anomalous high-latitude SST warming between 2.9 and 2.7 Ma, we simulated the injection  4) cross-section of the North Atlantic, respectively. We find that the MOW intrusion has a two-fold effect on North Atlantic hydrography: firstly, the zonal section ( Fig. 4a) indicates that the MOW plume settles dominantly at intermediate water depth (800-000 m) during its westward propagation leading to a warming of~2°C in the mid-Atlantic; secondly, the model results highlight the formation of a subsurface heat channel (Fig. 4b) that seemingly traverses the accumulated subtropical heat from intermediate water depths towards the surface of the high-latitude North Atlantic (>50°N) where a surface warming of~2-3°C can be observed-a process This subsurface heat transport strongly resembles results from modern drifter studies on the heat connectivity between subtropical and subpolar gyres. These results indicate that the communication between both gyres is at the surface level strongly limited due to a strong vorticity gradient 23,42 . As the vorticity gradient between both gyre systems decreases with increasing water depth the heat stored within the subtropical gyre can be channelled north along isopycnals that outcrop in the subpolar gyre 23 . The injection of heat by the MOW below the subtropical gyre deepens the subtropical thermocline and thickens the density layers through which subtropical waters reach the subpolar gyre. Thick isopycnals increase the subtropical-subpolar flow at levels not constrained by vorticity gradients (Fig. 2). Therefore, the MOW provides a channel for subtropical waters to reach higher latitudes and its effect can be assessed by the thermal gradient between the surface and deep thermocline at mid-latitudes. Modern hydrographic data show that the injection of dense MOW decreases the stratification of the mid-latitude North Atlantic and thus opens this subsurface heat channel, which supplies warmth to high latitudes (Fig. 2 41 North Atlantic (Fig. 2). Indeed, within the iMOW interval, at~2.8 Ma we find a warming of~3°C of intermediate water depths (~800 m) of the central North Atlantic IODP Site U1313 (Fig. 3f) 43 . This temperature increase argues for a heat accumulation at subsurface levels and thus strongly mimics the model predictions for the zonal MOW behaviour, also in terms of magnitude. In contrast, we observe a cooling of~4°C of intermediate waters (~1200 m water depth) at Site DSDP Site 548 during the iMOW interval, reaching lowest temperatures at~2.8 Ma (Figs. 2 and 3e). Under modern conditions, Site DSDP 548 is bathed in the northward propagating branch of the MOW 41 . Thus, the temperature decline argues for the absence or substantially reduced influence of MOW at this site during the iMOW interval. Based on our model results, this lack of MOW influence at Site DSDP Site 548 during the iMOW interval originates from a shoaling of the MOW plume, which is also in line with previous suggestions 41 . Hence, the modelled heat anomalies across the intermediate-depth North Atlantic generated by strong MOW production are in good agreement with proxy evidence. This is particularly noticeable when considering the limitations of the model to resolve spatially complex hydrographic conditions such as prevalent along the European margin, which is related to limited spatial resolution towards enabling long-term integrations for paleoclimatic applications. The ability to represent the impact of MOW on the oceanic heat distribution in a low-resolution model setup hints that the subsurface heat channel across the Atlantic Ocean may be a robust feature of past, present and future ocean circulation.
A long-term perspective of MOW influence on the North Atlantic. To assess the potential of this subsurface heat channel during other time periods, we analysed the existing mid-to highlatitude SST records back to 3.6 Ma when MOW started to intrude into the North Atlantic on a similar scale as today 34,41 . In fact, we find that when considering long-term trends, a similar anomalous high-latitude warming phase in the North Atlantic between 3.5 and 3.1 Ma dovetails strongly increased salinities in the western Mediterranean Sea (Fig. 5) as well as coincides with similar timed erosional features along the continental slope of the Gulf of Cadiz 34 . This provides strong indications that an enhanced MOW production and development of a subsurface heat channel between the mid and high latitudes in the North Atlantic might have been active during the glacial period M2 (~3.3 Ma) as well as during the high pCO 2 period KM5c (3.21-3.19 Ma) 44 . However, a detailed study of this interval, in contrast to the 2.9-2.7 Ma interval, is currently not possible due to poor preservation of the sedimentary sequence from the Gulf of Cadiz prior to~3.2 Ma as well as the lack of intermediatewater-mass records from the mid-latitude North Atlantic.
An active subsurface heat channel driven by strong MOW production during KM5c, considered an analogue for future climate change 45 , might nevertheless provide a blueprint of its impact on the North Atlantic hydrography under increasing pCO 2 levels. Model results show that during KM5c the Mediterranean experienced 2-3°C warmer temperatures, shifting its hydrological balance into greater net loss of freshwater 46 . This is a similar climate development as predicted for the Mediterranean region until the year 2100 47,48 . Consequently, the lack of freshwater influx to the Mediterranean during KM5c probably boosted its salinities and with it MOW production (Fig. 5). The increased MOW production then might also have led to a rerouting of heat to the high-latitude North Atlantic via the subsurface aligning with established SST dynamics for this time period 8 . Our mechanism thus provides an explanation for the observed high-latitude SST anomaly during KM5c, which has been commonly attributed to polar amplification instead 45 .
In summary, we argue that a MOW-driven subsurface heat channel might have played a role in delaying Northern Hemisphere ice sheet growth against the backdrop of global cooling and declining pCO 2 whilst also holding insights into the SST evolution of the North Atlantic under future climate scenarios. As aridity is projected to increase severely in the circum-Mediterranean region until the year 2100 47 , we expect MOW production to increase in the future. This in turn would invigorate the subsurface heat transport to the high-latitude North Atlantic, as indicated for the interval of KM5c, and thus potentially further amplify the effects of global warming in this climatically extremely sensitive region.   20,49 . The LIW is formed in the eastern Mediterranean Sea in the highly evaporative Levantine Basin during wintertime cooling. It feeds into the uWMDW that originates in the Gulf of Lyons and constitutes the densest portion of the MOW 20 . During strong North African aridity phases, the fresh water supply to the Mediterranean, e.g. by strong Nile discharge, is supressed, increasing the salinity of the Mediterranean intermediate-and deepwater masses and in turn enhancing MOW formation 28,30,31 .

Methods
To reconstruct the salinity of the uWMDW as a measure of MOW production, we utilized the existing stable oxygen isotope (δ 18 O) and Mg/Ca-based bottom water temperatures (BWT) derived from benthic foraminifera Cibicides spp. from Site ODP 978 in the Alboran Sea 41 . Site ODP 978 is located in the Alboran Sea north of the Al-Mansour Seamount (36°13ʹN, 02°03ʹW; Fig. 2) at~1920 m water depth, and is thus bathed in uWMDW 41  Modelling the influence of MOW on the North Atlantic water column. The COSMOS PlioMIP2 core simulation Eoi400 48 was chosen as a starting point for studying the impact of MOW on the North Atlantic hydrography during the Plio-Pleistocene transition 48 . It is based on the state-of-the-art reconstruction of mid-Piacenzian paleoenvironment by 53 . Consequently, it provides a reference setup that considers best knowledge on the mid-Pliocene state of ice sheets, configuration of continents and ocean gateways, as well as carbon dioxide concentrations. However, the COSMOS PlioMIP2 core simulation does not include MOW injection into the North Atlantic and, furthermore, the topographic details of the Strait of Gibraltar are limited due to the coarse spatial resolution of the ocean model (~3.0°× 1.9°formally at a bipolar curvilinear model grid). Hence, to simulate the intrusion of MOW and to analyse its impact on North Atlantic circulation and hydrography, we implemented a prolonged Newtonian salinity restoring towards a reference salinity of 45 PSU at the Strait of Gibraltar (34°N, 352°E) at intermediate water depth (about 700 m). The 45 PSU input MOW salinity is derived from converting the maximum δ 18 O sw-ivf from ODP Site 978 (Alboran Sea) into bottom-water salinity (PSU) by applying the modern Mediterranean Sea-specific δ 18 O sw-ivf -salinity relationship following 54 . We argue that because Late Pliocene background conditions are comparable to modern in terms of paleogeography and climate 41,55 , the applied modern Mediterranean Sea δ 18 O swivf -salinity relationship can be used as an approximation. This assumption is supported by numerical model simulations for the Late Pliocene that suggest a similar δ 18 O sw-ivf -salinity relationship between Modern and Late Pliocene, and even highlight an increase in Mediterranean surface salinity under Late Pliocene  41 e Sapropels (representing wet phases) in Eastern Mediterranean Site 969 67 . The grey shadings indicate intervals of anomalous SST warming in the high-latitude North Atlantic concomitant to inferred increases in MOW production, including the "iMOW-Phase" discussed in the text.
boundary conditions as proposed in this study 55 . Notably, the modelled global ocean δ 18 O sw-ivf -salinity relationship for the Late Pliocene has an even smaller slope (0.22) 55 than the modern one (0.27) 54 , which would generate even higher salinities for the MOW in comparison to our conservative estimates based on the modern δ 18 O sw-ivf -salinity relationship 54 . Based on our model simulations, the denser MOW flow core during the iMOW interval settles at roughly~1200 m water depth within the North Atlantic (Fig. 4), which is deeper than its modern main flow core depth between 500 and 800 m water along the continental slope of the Gulf of Cadiz 28,30 . Notably, the simulated model depth is in line with MOW flow core settling depth reconstructions from the Last Glacial Maximum when a MOW with higher-than-modern salinity (<40 PSU) settled at water depths of roughly 800-1200 m 56 . Hence, the simulated heat transport already does take into account the deeper MOW flow path as a consequence of the increased salinity.
Starting point of this modelling study is a quasi-equilibrium mid-Pliocene climate state derived from the COSMOS PlioMIP2 mid-Pliocene core simulation Eoi400 48 . We initialized the model from the mid-Pliocene core climate state after a model spin-up over 1400 model years. The climate simulation was integrated thereafter with Newtonian salinity restoring as described above. The strength of the AMOC reacted to the implied more saline MOW with an overshoot after about 1750 model years. The impact of the mimicked MOW salinification on the hydrography of the North Atlantic Ocean is investigated based on this AMOC overshoot.
COSMOS model description. The Community Earth System Models (COSMOS) provide a numerical representation of the coupled atmosphere-ocean-land system with explicit consideration of the land-surface carbon cycle, including dynamic vegetation. The climate system is represented by the atmospheric general circulation model ECHAM5 57 with river runoff scheme, the land surface and carbon cycle model JSBACH 58 with vegetation dynamics 59 , and the ocean general circulation model MPIOM 60 with a dynamic-thermodynamic sea ice scheme based on the model by 61 . The COSMOS has been applied with a resolution of T31 (3.75°× 3.75°) and 19 hybrid sigma-pressure levels in the atmosphere and a bipolar curvilinear GR30 grid with formal resolution of~1.8°× 3.0°and 40 z-levels in the ocean domain. Ocean and atmosphere are coupled once per model day via the OASIS3 coupler 62 . Time steps of atmosphere and ocean are 2400 and 8640 s, respectively. Methods employed in COSMOS towards improvement of ocean circulation at coarse model resolution are outlined in the literature; for an overview see 63 as well as 48 , who present a detailed description of the model and of its application in the framework of paleoclimate research. The model does not necessitate any flux corrections. Beyond the implementation of Newtonian salinity restoring at the Strait of Gibraltar towards mimicking more saline MOW, we have not applied any artificial fluxes in the model setup.
Simulating the impact of MOW in the COSMOS model. Towards mimicking more saline MOW in the model, we perform a localized restoring of salinity (S) towards a target value S ref that is chosen to represent the characteristics of MOW (see above). The parameterization of high-saline Mediterranean Outflow via Newtonian salinity restoring is performed at one grid cell at a depth of about 700 m at the Strait of Gibraltar (34°N, 352°E) following Eq. 1.
The ocean model locally restores salinity towards the, in comparison to the simulated salinity, increased target value (S ref = 45.0) by computing a new salinity (S new ) that is based on the salinity S old of the previous time step and a restoring time-constant (f) of 10.368 × 10 6 s −1 that refers to a time period of 4 months.
Overshoots in the AMOC, which result from the application of Newtonian salinity restoring, represent a surrogate to the impact of more saline MOW on large-scale circulation and oceanographic patterns of the Pliocene Atlantic Ocean.

Data availability
All referenced data sets are online available via pangea.de. The model data of the COSMOS PlioMIP2 mid-Pliocene core simulation and the respective pre-industrial control state are available from a repository at the University of Leeds (ftp://see-gw-01. leeds.ac.uk).