Climate‐Induced Variability in Mediterranean Outflow to the North Atlantic Ocean During the Late Pleistocene

Mediterranean Outflow Water (MOW) adds salt and density to open ocean intermediate waters and is therefore an important motor of Atlantic meridional overturning circulation (AMOC) and climate variability. However, the variability in strength and depth of MOW on geological timescales is poorly documented. Here we present new detailed records, with excellent age control, of MOW variability from 416 ka to present from rapidly accumulated marine sediments recovered from the West Iberian Margin during Integrated Ocean Drilling Program (IODP) Expedition 339. Our records of X‐ray fluorescence (XRF), physical grain size, and paleocurrent information from the anisotropy of magnetic susceptibility (AMS) indicate (i) a close relationship between the orientation of principle AMS axes and glacial‐interglacial cycles and (ii) two distinct regimes of MOW behavior over the last ~416 kyr in grain‐size and AMS variability at orbital (mainly precessional) and suborbital timescales. Between marine isotope stage (MIS) 10 and MIS 4, MOW was focused at a generally shallow depth on the West Iberian Margin, and changes in MOW strength were strongly paced by precession. A transition interval occurred during MISs 5 and 4, when MOW deepened and millennial‐scale variability in flow strength was superimposed on orbitally paced change. During MIS 11 and from MIS 3 to present, MOW was deeply focused and millennial‐scale variability dominated. We infer that late Pleistocene variability in MOW strength and depth were strongly climate influenced and that changes in circum‐Mediterranean rainfall climate were likely a primary control.


Introduction
The relationship between Mediterranean Outflow Water (MOW) and North Atlantic circulation strength and structure is complex. MOW provides a salt source for intermediate depths of the North Atlantic, but the effects of this input are not well understood. It is widely (but not universally) suggested that the salt supplied by MOW significantly increases the rate of North Atlantic overturning today and had variable impacts during the last glacial cycle Llave et al., 2006;Rogerson et al., 2010Rogerson et al., , 2012Swingedouw et al., 2019;Voelker et al., 2006). MOW probably exerts significant circulatory impact following Greenland stadials, forming a crucial part of the negative feedback mechanism that promotes the reinstatement of vigorous North Atlantic deepwater formation and overturning from its weakened cold-interval state (Llave et al., 2006;Voelker et al., 2006). Furthermore, significant reorganizations of shallow and intermediate water mass structure in the North Atlantic may arise in response to changes in MOW-injection depth  and the strength of the Azores Current, and consequently, European ice sheet dynamics are sensitive to MOW entrainment into North Atlantic intermediate water masses . In contrast to the idea that increased MOW-fed salinity to the North Atlantic enhances overturning circulation, one recent numerical analysis of the impact of sapropel events suggests that the upper AMOC cell may weaken or strengthen, depending on the rate of freshwater input to the Mediterranean (Swingedouw et al., 2019).
On millennial to orbital timescales, MOW variability is primarily controlled by the interplay between the dimensions of the Strait of Gibraltar, the density contrast between Atlantic and Mediterranean water masses, and the density structure within the North Atlantic, and these factors affect both the strength of outflow and isopycnal equilibration (Rogerson et al., 2005(Rogerson et al., , 2010(Rogerson et al., , 2012. Glacioeustatic sea-level change was the primary control on MOW over the last glacial cycle (Rogerson et al., 2012). The~120 m sea-level fall (relative to present) associated with the Last Glacial Maximum severely reduced the size of the channel at the Strait of Gibraltar, decreasing exchange between the Mediterranean and Atlantic. The result was a less voluminous but denser and deeper-equilibrating MOW relative to present (Rogerson et al., 2012). Changing density gradients across the Strait of Gibraltar exerted a secondary control on MOW flow strength over the last glacial cycle, modulated by precipitation-evaporation around the circum-Mediterranean region, riverine input, and inflow of Atlantic water along the surface at the Strait of Gibraltar (Rogerson et al., 2012). Freshwater input into the Mediterranean was less pronounced during stadials and precession maxima (insolation minima, i.e., aphelion near the northern summer solstice) than during interstadials and precession minima (insolation maxima, i.e., perihelion near northern summer solstice) because of diminished northward migration of the Intertropical Convergence Zone (ITCZ) and African summer monsoons, and, in winter, steering of storm tracks (Bosmans et al., 2015;Kutzbach et al., 2014) or changes in sea-level pressure at the Azores High (Bosmans et al., 2020). The resultant combination of more saline surface waters and cooler air temperatures acted to enhance Levantine Intermediate Water (LIW) and deepwater formation in the eastern Mediterranean and, in turn, the rate of Mediterranean Outflow and salinity export to the North Atlantic (Bahr et al., 2015;Rogerson et al., 2012;Toucanne et al., 2012). This mechanism is hypothesized to have helped to reinstate vigorous overturning in the North Atlantic following stadial surface water freshening (Voelker et al., 2006).
Thanks to Integrated Ocean Drilling Program (IODP) Expedition 339 (Mediterranean Outflow) (Expedition 339 Scientists, 2012), several studies have explored MOW evolution deeper into the Pleistocene. Benthic foraminiferal assemblages over the last 900 kyr show distinct glacial-interglacial contrasts at Site U1391, and these have been related to changes in MOW activity (interglacial maxima and glacial minima with peak velocities during marine isotope stage [MIS] 11) (Guo et al., 2017). These findings are consistent with the suggestion that higher salinity MOW and a reduced Atlantic vertical density gradient during glacials drive deeper isopycnal equilibration of MOW and therefore a reduced contribution to Atlantic overturning (Rogerson et al., 2005). At IODP Site U1386 in the Gulf of Cadiz, δ 18 O and grain size (mean weight percent of the 150-63 μm fraction) reveal three distinct phases of MOW variability back to 570 ka (Kaboth et al., 2017). Significant MOW variability was identified between 570 and 475 ka (Phase III) and from~130 ka to the present (Phase I). The grain-size record of Site U1386 suggests that MOW was significantly slower in the intervening period between 475 and 130 ka (Phase II). Shifts in MOW water mass sourcing are invoked to explain these changes, but records are needed from sites in different water depths to test this hypothesis, to determine the history of MOW depth equilibration variability and to understand the forcing mechanisms involved.
Here we present new multiproxy records from a two-site depth transect on the West Iberian Margin to study both the strength and depth variability of MOW on millennial to orbital timescales over the past 416 ka and to make inferences about the relative importance of its driving mechanisms through time. We find that MOW behavior shows a strong imprint of climate-induced variability with two distinct regimes of orbital and millennial timescale variability.

Regional Modern Hydrography
On short timescales (years), the composition of MOW is highly variable. Any combination of four water masses, LIW, Tyrrhenian Deep Water (TDW), Winter Intermediate Water (WIW), and Western Mediterranean Deep Water (WMDW), can contribute to changes in MOW composition (Millot, 2009). However, it is generally agreed that over glacial-interglacial timescales, LIW is the most prominent contributor (~90%) to compositional change. LIW forms in the eastern Mediterranean when surface water temperatures are lowered and salinities increase in response to strong outbreaks of cold dry air masses over the region during winter . LIW then flows westward to the Strait of Gibraltar. Across the Alboran Sea the WMDW shoals from~800 up to 284 meters below sea level (mbsl) because of Bernoulli aspiration with the fast flowing overlying LIW promoting increased buoyancy in WMDW, causing it to be drawn upwards over the Camarinal Sill (the shallowest point in the Strait of Gibraltar). Upon exiting the Strait of Gibraltar, the relatively cold and saline MOW flows downslope at 0.57 ± 0.26 Sv and splits into an upper limb and a lower limb (Millot, 2009;Tsimplis & Bryden, 2000). These limbs are formed in response to gradual entrainment of Atlantic water into the outflowing MOW plume. The outer plume undergoes frictional mixing and slows relative to the central core (O'Neill-Baringer & Price, 1999). Ekman veering causes the central unmixed flow to veer downward at a steeper angle and to be injected into the Atlantic at a greater depth than the surrounding outer plume Mediterranean water (Rogerson et al., 2012). While the two limbs have similar salinities,~38.5 psu, there is a 1.5-3.5°C temperature difference between them causing them to attain different isopycnal depths between~600 and~1,500 m. In the Gulf of Cadiz, MOW can be traced flowing to both the west and south; however, a significant proportion turns northwards, flowing around the West Iberian Margin (Iorga & Lozier, 1999) (Figure 1). Here, further from the Strait of Gibraltar, greater entrainment of Atlantic water means that the two limbs are less distinct than in the Gulf of Cadiz. Thereafter, MOW can be traced as far north as the Norwegian-Greenland Sea (Lozier & Stewart, 2008).
On the West Iberian Margin, the seasonally variable Portugal Current system occupies the surface down tõ 100 m water depth. Seasonal shifts in strength of the Azores High pressure system promote direction-offlow change with the northward flowing Iberian Poleward Current dominant during winter and the southward flowing Portugal Current dominant during summer (Relvas et al., 2007). Below this surface system there are two components of the Eastern North Atlantic Central Water (ENACW): a northward flowing subtropical component (ENACWst) at~200-300 m and a southward flowing colder nutrient-rich subpolar component (ENACWsp) at~300-400 m (Arístegui et al., 2006). MOW flows northward beneath ENACWsp from ~600 down to~1,500 m. North Atlantic Deep Water (NADW) dominates below MOW and flows southward. Antarctic Bottom Water (AABW) occupies depths below the NADW, but can increase in volume, shoaling to reach~2,500 mbsl during glacials (Hodell et al., 2013).

IODP Sites U1391 and U1385
Contourites are deep-sea sediments deposited under the influence of a bottom water current. MOW-generated contourite sediment packages extend from the Gulf of Cadiz to the West Iberian Margin. IODP Site U1391 is the most distal MOW-entrained site from the Strait of Gibraltar drilled during Expedition 339 (Figure 1). The stratigraphic interval studied here consists of sandy contourites and alternating red/brown and green/gray calcareous muds (Expedition 339 Scientists, 2013b). Sedimentation rates at Site U1391 are an order of magnitude higher than at typical pelagic sites, averaging~27 cm/kyr over the last 1.5 Ma (Hernandez-Molina et al., 2014). While there is no fixed boundary between the upper and lower limbs of MOW on the western Iberian Margin, U1391 is situated closer to the lower limb in terms of its water depth (1,085 mbsl).
IODP Site U1385 is situated northwest of Site U1391 ( Figure 1) in deeper water (~2,585 mbsl), under the influence of NADW and southern sourced waters rather than MOW. Sediments at Site U1385 comprise nannofossil muds and clays with cyclical carbonate and color variability (Expedition 339 Scientists, 2013a) and average sedimentation rates of~11 cm/kyr . Site U1385 was drilled near the position of Core MD01-2444, which yielded planktic and benthic oxygen isotope signatures closely resembling temperature records from Greenland and Antarctic ice cores, respectively (Hodell et al., 2013;Shackleton et al., 2004). Marine sediments from the Iberian Margin also provide detailed records of European terrestrial climate because of the narrow continental shelf and proximity of the Tagus River (Tzedakis et al., 2009(Tzedakis et al., , 2015. The clear and unequivocal links that Site U1385 stratigraphy provide between the marine, terrestrial, and ice core records make it an ideal reference section at which to develop a high-resolution chronology .

Samples
A total of 106 u-channel samples (with typical dimensions of~150 × 2 × 2 cm 3 ) was taken from the archive-half core sections, following the Site U1391 shipboard stratigraphic splice. The u-channel samples cover a continuous sediment sequence of~120 m with some overlapping intervals between samples from different holes. In addition, a total of 274 paleomagnetic cubes (each with a volume of~7 cm 3 ) and 37 bulk sediment samples was collected at selected depth levels from the same~120 m sediment sequence.

X-ray Fluorescence Scanning
Site U1391 u-channels were scanned on the ITRAX X-ray fluorescence (XRF) core scanner at the British Ocean Sediment Core Research Facility (BOSCORF) at the National Oceanography Centre Southampton. A surface profile scan was taken to allow detector height adjustments downcore to ensure constant detector to sediment distance. A 4-μm-thin SPEXCertiPrep Ultralene film covered the sediment surface before scanning. Data were collected at 0.5 cm intervals with a 15 s count time at 30 mA current and 30 keV voltage. Split archive-half core sections from Holes U1385D and U1385E were scanned at 1 cm resolution on an Avaatech core scanner at the Godwin Laboratory at University of Cambridge, as reported by Hodell et al. (2015).
The ratio of calcium to titanium (Ca/Ti) is commonly used as a correlative tool between cores in marine sediments (and is used as such in this study) and reflects the changing relative proportions of biogenic and detrital sediment sources at Site U1385 (Hodell et al., 2013). Zirconium is commonly found in resistate minerals, while rubidium is more highly concentrated in clay minerals; therefore, the ratio between the two (Zr/Rb) is used as a grain-size proxy Rothwell et al., 2006). Element counts in XRF data not only depend on concentration but also on physical properties, such as water content, geometry of the sample, and programmed settings of the scanner. Using log ratios has been shown to be an effective method of minimizing these effects, allowing elemental variations to be compared between records (Weltje & Tjallingii, 2008).

Anisotropy of Magnetic Susceptibility Measurements
The fabric of paramagnetic and ferrimagnetic grains in a sample can be estimated by measurement of the anisotropy of magnetic susceptibility (AMS). AMS describes the directional variation of magnetic susceptibility and is expressed as a second-rank tensor converted to a magnitude ellipsoid (King & Rees, 1962). Three eigenvalues are used to represent the maximum, intermediate, and minimum axes of the ellipsoid, (Hrouda, 1982). The inclination and declination of each principal axis of the ellipsoid can be plotted on lower hemisphere projections and, alongside other calculated parameters, can reveal information about the nature of deposition ( Figure 2). In the majority of cases where a current is acting on sediment, the orientation of the K max axis is parallel to the flow direction (Taira, 1989) (Figure 2c).
Other key AMS parameters are summarized in supporting information Table S1.
AMS of the U1391 paleomagnetic cubes was measured on a Kappabridge KLY-4S in the paleomagnetism laboratory at the University of Southampton. Susceptibility was measured in 64 different directions while the sample was slowly spinning. Declinations of each AMS axis were corrected using the vector-mean paleomagnetic declination of each core. This was calculated from natural remnant magnetization (NRM) measurements of the u-channels from the archive-half core sections, taking the mean of 12 AF Gravitational settling of grains in quiet water settings results in well-spread K max axes declinations and K min clustered around the center, as shown in (b2) stereographic projection. (c1) Intermediate strength current flow imbricates K max axes opposite to the flow direction and skews K min axes from vertical into the direction of flow. In this instance K max axes are grouped to the north, and K min axes skewed to the south, both of which indicate southward flow, as shown in (c2). (b2-c2) Lower hemisphere projections of AMS K max , K int , and K min declination and inclination values from two intervals demonstrating each respective fabric type. The figure is based on similar figures in Taira (1989) and Tauxe (2010). demagnetization steps from 20 up to 90 mT (Xuan, pers. comm. 2019). U1391 cores relevant to this study are typically 9.5 m long covering a total of~30 kyr (assuming a mean sedimentation rate of 27 cm/kyr), over which the geomagnetic field should approximate a geocentric axial dipole field with secular variations averaged out, giving an expected declination of 0°for Brunhes-aged sediments at the site. Therefore, the vector-mean NRM declination of each core provides an estimate of the rotation needed to restore the azimuthal orientation of the core relative to north. NRM declinations of some cores show an apparent along-core trend, indicating that the core has twisted during the coring process (e.g., Lanci et al., 2004). For cores with an apparent trend in NRM declinations, we estimate a linear trend line and calculate declination correction values for cube samples from the core based on their relative position within the core.

Physical Grain-Size Analysis
Physical grain size is well established as an indicator of paleocurrent speed variations (McCave et al., 1995(McCave et al., , 2017. The size fraction 10-63 μm represents the sediment that is hydrodynamically sorted and therefore is most representative of current strength. The mean of this range is known as the sortable silt mean proxy (SS) and has been linked to flow speed (McCave et al., 1995(McCave et al., , 2017. Thirty-seven bulk samples from Site U1391 were weighed, treated with 0.2% sodium hexametaphosphate solution, left on a shaker table overnight to disaggregate, and then sieved at 106 μm. The samples were then washed with 10% HCl to remove carbonate minerals and calcareous microfossils. Physical grain-size sample processing often includes a biogenic silica removal step. Inspection of smear slides across a range of climatic intervals showed that biogenic silica is extremely scarce (<1%) in Site U1391 sediments, consistent with shipboard observations (Expedition 339 Scientists, 2013b). The biogenic silica removal step was therefore not included during sample processing. Physical grain-size distributions of the processed samples were measured using a Multisizer 3 Coulter Counter at the University of Cambridge Godwin Laboratory. Samples were suspended in an electrolyte solution, treated with 2 minutes of ultrasound, and stirred. The suspended material was then passed through an aperture of 200 μm where volume-related impedance changes are caused by displacement of the electrolyte as the grains pass through. Grain-size distribution was measured in log 2 steps in the range between 5 and 100 μm for a total of 20,000 counts. Measurements for each sample were repeated at least three times and typically six times to increase statistical reliability.

Chronology and Sedimentation Rate
Our new XRF core-scan records from Site U1391 are presented and compared to data from Site U1385, Site MD01-2444, and the synthetic Greenland oxygen isotope record in Figure 3. Over most of the shared stratigraphy, the XRF records from Sites U1391 and U1385 are remarkably similar to one another, allowing detailed correlation, in particular between ln(Ca/Ti). In general, ln(Ca/Ti) ratios from the two sites also show close similarity to an alkenone-based record of sea surface temperature from Site MD01-2444 on the Iberian Margin (Martrat et al., 2007). The chronology for Site U1391 was constructed by transferring the well-established age model from nearby Site U1385 , using ln(Ca/Ti) as the correlation tool ( Figure 3 and Table S2). Correlation of ln(Ca/Ti) between the two sites is mostly straightforward, except in the intervals 65-73, 140-190, and 348-380 ka, where the correlation is less clear. Planktic and benthic δ 18 O records from Site U1391 (Guo et al., 2017) were placed on our constructed age model and compared with equivalent U1385 records ; Figure S1). This figure confirms the robustness of our ln(Ca/Ti) correlation between Sites U1391 and U1385 because major climate transitions and much of the millennial scale variability documented in the δ 18 O records are well aligned between the two sites.
The synthetic Greenland age model from Site U1385 was transferred to Site U1391 and is based upon the correlation of ln(Ca/Ti) at Site U1385 to the same parameters obtained from nearby piston Cores MD01-2444 and MD01-2443 . Hodell et al. (2013) correlated sea surface temperatures and planktic δ 18 O records at MD01-2444 and MD01-2443 to the synthetic Greenland temperature record (Barker et al., 2011) (see Figure 3). The synthetic Greenland temperature record was constructed using a bipolar seesaw model prediction, making use of the relationship between Antarctic and Greenland temperature variations, and was chronologically tied to the absolutely dated Chinese speleothem record (Barker et al., 2011).
From 416 to 130 ka at Site U1391, the sedimentation rates at Site U1391 remain steady, at between~28 and 30 cm/kyr (Figure 4). Between~130 and~50 ka the implied sedimentation rates increase steadily up to 39 cm/kyr before decreasing back to around 24 cm/kyr in the Holocene.

AMS as a MOW Indicator at Site U1391
AMS parameters from Site U1391 cube samples, including K max declination, degree of anisotropy (P j ), and q values (see Table S1 for AMS parameters explanations), are shown in Figure 5. K max declination shows distinct behavior throughout. Prior to~130 ka, K max axes orientations show apparent glacial-interglacial variation, with the majority of the orientations grouped around the southeast during glacials, particularly between~320 and 130 ka ( Figure S2). After~130 ka, changes in K max declination became dominated by millennial-scale variability.
α 95 values associated with the K max declination data are mostly less than 10°and are displayed in a histogram in Figure S3. A lower hemisphere projection of K max , K int , and K min declination and inclination values from all samples is presented in Figure S4. Throughout the record K min inclinations are predominantly clustered perpendicular to the horizontal with a mean value of 88.9°( Figure S4). K max declinations, across the entire studied interval, are not only grouped to the southeast and northwest but also spread in directions between these two groups. P j -T and q-β plots ( Figures S5 and S6), however, indicate that the fabric is mostly oblate and the sediments are likely current influenced. These P j -T and q-β values indicating current activity are not restricted to certain K max declination directions. Influences such as a slope at the depositional site, coring, or sampling disturbance could have led to deviations in K max declination. It is also possible that downslope movement gave rise to more randomized directions in certain intervals, although this possibility is less likely according to the P j -T and q-β plots ( Figures S5 and S6). Furthermore, a stronger current may lead to grains deposited with a K max declination perpendicular to flow, though this is suggested to mainly occur in non-Newtonian fluids (Novak et al., 2014;Tauxe et al., 1998) and unlikely to be important at Site U1391 because the diagnostic behavior of "streak-distributed" K min axes in the direction of flow on lower hemisphere projections (Park et al., 2013) is not observed.
Strong inverse correlation is observed between the degree of anisotropy (P j ) and ln(Ca/Ti) ( Figure 5). Orbital scale variation in ln(Ca/Ti) is linked to dilution of carbonate by terrigenous input from the shelf on the West Iberian Margin, while productivity changes largely reflect millennial scale variability (Hodell et al., 2013). Changes in P j have been linked to both flow speed and delivered sediment composition (Yokokawa & Franz, 2002). Figure 5 shows that when P j is low, the magnetic fabric, q, and physical grain-size proxies indicate stronger MOW. Given that ln(Ca/Ti) is high in these intervals, generally corresponding to a warmer climate (Figure 3), it is also likely there is a greater degree of bioturbation that may lower P j . When P j is high, lower ln(Ca/Ti) indicates a greater proportion of terrestrial input. Directional AMS data do not appear to have any relationship with the degree of anisotropy. We infer, therefore, that P j at Site U1391 mainly represents the foliation of the sediment produced at times of higher terrestrial input, lower bioturbation, and slower flow, rather than a current-induced lineation.
From MISs 10 to 4, K max declination variability shows remarkable first-order correlation with relative sea level ( Figure 5). During glacial interval MISs 6, 8, and 10, the majority of the K max declination values dip to the southeast, suggesting a predominantly northwest imbrication and therefore northwestward flow that is likely attributable to MOW (see Figures 2, 5, and 6). During most of interglacial MISs 5, 7e, 9a, 9e, and 11, however, K max declination is more scattered, with many samples in these intervals showing clear dips to the north suggesting a degree of southward imbrication (Figures 2, 5, and 6). This observation implies a less stable and possibly weaker southward flowing water mass over Site U1391 at these times than during glacials. Given changes in Atlantic vertical density structure and a less dense MOW during interglacials (Rogerson et al., 2012), MOW likely shoaled during these interglacial stages and may only have been active in water depths shallower than Site U1391. The water mass lying underneath MOW is the southward flowing, eastern NADW (ENADW). If the MOW shoaled during these interglacials it follows that ENADW also shoaled and/or expanded to influence sedimentation at Site U1391 at these times.

Orbitally Paced MOW Influence at Site U1391 During MISs 10 to 4
Our grain-size data show predominantly polymodal and moderately sorted distributions. The geometric mean grain size at Site U1391 is a coarse silt of~20.8 μm, with the highest value of~30.5 μm occurring at 139 ka. Although these data are of relatively low temporal resolution, the mean values of grain size over time have similar patterns of variability to those documented in our high-resolution ln(Zr/Rb) data series ( Figure 5). The positive correlation between mean physical grain size and ln(Zr/Rb) at Site U1391 ( Figure S7) suggests that ln(Zr/Rb) functions well here as a mean physical grain-size proxy, as discussed in section 3.5 and suggested elsewhere Rothwell et al., 2006). At Sites U1391 and U1385 ln(Zr/Rb) shows precession-paced variability between MISs 10 and 4, while in the absence of high-amplitude precession forcing, millennial-scale variabilities dominate in MIS 11 and MISs 3 to 1 (Figures 5 and 7). After MIS 5, millennial-scale cold events (vertical gray bars in Figure 9) are marked by distinctly coarser sediments at Site U1385 compared to Site U1391.
Our grain-size records from Sites U1391 and U1385 rule out a major influence by sea level-forced downslope movement. R 2 coefficients between Site U1391 ln(Zr/Rb) and downcore variation of each grain-size bin population are shown in Figure S8. The ln(Zr/Rb) record correlates best with grain sizes~30-70 μm, mostly in the sortable silt range, which is therefore more likely hydrodynamically sorted (MOW) material. The sortable silt proxy has been used previously on the Iberian Margin to identify paleocurrent changes (Hall & McCave, 2000). Finer sediments (grain sizes <~20 μm), on the other hand, which are more likely to be controlled by lower energy downslope transport, appear to be inversely correlated with the ln(Zr/Rb) data ( Figure S8). Further evidence of current speed-controlled grain-size variability is found in the positive correlation between grain size and q values (magnetic fabric, see Table S1). q values tend to peak at precession maxima ( Figure 5), indicating greater importance of shear stresses acting tangentially along the bed (i.e., a stronger current influencing deposition) at these times.
The clear precession signal in the ln(Zr/Rb) grain-size proxy (Figures 5 and 7) from MISs 10 to 4 likely reflects changes in flow velocity at Site U1391 and points to the influence of circum-Mediterranean rainfall climate variability. Reductions in both winter and summer precipitation at insolation minima (Bosmans et al., 2015(Bosmans et al., , 2020Kutzbach et al., 2014;Trauth et al., 2009) likely contributed to increased LIW formation and outflow strength. Our ln(Zr/Rb) proxy record of current speed shows remarkable alignment with changes in spring insolation at 35°N (Figure 5). This result suggests a stronger response of MOW to late winter/spring climate forcing of eastern Mediterranean hydrography at these timescales than to summer insolation.
Although Site U1385 lies below the expected depth of MOW influence, it records a pattern of ln(Zr/Rb) variability similar to that observed at Site U1391 (Figure 7). Additionally, the insolation forced grain-size pattern continues through MISs 5, 7, and 9, when the AMS K max declination data indicate that MOW had shoaled above Site U1391 ( Figure 5). We infer that terrestrial sediments deposited on the West Iberian Margin (Hodell et al., 2013;Thomson et al., 1999) were hydrodynamically sorted by MOW in transit to sites below MOW (Figure 8). Similar hydrodynamic sorting effects on organic carbon transport have been inferred for a sediment core near the location of Site U1385 (Magill et al., 2018).
During the precession maximum between~160 and 170 ka, mean physical grain size and ln(Zr/Rb) at Site U1391 do not coarsen as they do in previous precession cycles, while ln(Zr/Rb) at Site U1385 does coarsen (highlighted in Figure 7). Following the grain-size hydrodynamic sorting mechanism (see Figure 8), we infer that MOW may have deepened below Site U1391 during this interval but still influenced grain size at Site U1385 through hydrodynamic sorting of terrestrial particles moving downslope to Site U1385. Notably, the K max declination does not vary greatly during this time interval (Figure 5), possibly indicating a partial presence of MOW that is sufficient to orient the magnetic grains but not enough to hydrodynamically sort the grain size. Other intervals of different ln(Zr/Rb) behavior, highlighted in Figure 7, tend to occur around cold events, suggesting that similar episodes of MOW deepening to depths below Site U1391 may have occurred many times during the last 416 ka.

Two Distinct Regimes of MOW Behavior on the West Iberian Margin
A first-order feature of our records is the strength of the orbitally paced (mainly precession) signals between MISs 10 and 4 (section 4.3). From MIS 5 we see strong millennial-scale variability, which is initially superimposed on the orbital signal, but from MIS 3 to present, the millennial-scale signal dominates (Figures 5  and 9). Thus, our records describe two regimes of temporal change and a transition interval between them. Regime I occurred from MIS 3 to present and during MIS 11. Regime II occurred between MIS 10 and MIS 5,  (Martrat et al., 2007;orange), ln(Zr/Rb) from Sites U1391 (blue) and U1385 (orange), and spring insolation at 35°N (Laskar, 2004;black). Cold events are marked by gray bars. Note that y-axis direction of the K max declination is reversed compared with Figure 4. inclusive. Through Regime II, MOW was orbitally paced and the average depth was shallower; it was present at Site U1391 during glacials and migrated above U1391 at peak interglacials. MOW behavior then underwent a transition during MISs 5 and 4, when the influence of both orbital and millennial-scale forcing is documented and the MOW plume deepened on the West Iberian Margin. Regime I began around 60 ka; MOW was on average deeper, migrating below Site U1391 during stadials, and responded to millennial-scale variability more strongly, as recorded by both grain size and K max declination signals. The timings of the position of MOW relative to Site U1391 in this reconstruction are summarized in Table S3.
The two distinct regimes interpreted here are consistent with the independently calculated sedimentation rate history of Site U1391 shown in Figure 4. When MOW is interpreted to either be present at or lie just above U1391 (Regime II), there is a steady sedimentation rate of~29 cm/kyr. Sedimentation rates increase during the transition interval consistent with MOW movement to a deeper position on the West Iberian Margin, with an average position closer to U1391, focusing more deposition at the site. The subsequent decrease in sedimentation rate around 50 ka is consistent with MOW being frequently located below and therefore bypassing the site.
For the last~130 ka, K max declination and sea surface temperature records show generally opposing behavior, declination implying northward flow at Site U1391 during warm intervals before flipping to indicate southward flow during stadials. This relationship contrasts with the one documented for orbital forced insolation minima and maxima during the preceding~235 kyr ( Figures 5 and 9 and section 4.2). As illustrated in the ln(Zr/Rb) records in Figure 9, during the transition interval from~130 to 60 ka, there is a strong imprint of millennial-scale events superimposed on the precession-paced grain-size cycles. From 60 ka, the precession signal disappears to produce a regime dominated by millennial change. During the latest~106 kyr, ln(Zr/Rb) records from Sites U1385 and U1391 diverge from one another during cold stadials (vertical bars in Figure 9) but resemble one another more closely during warmer interstadials. During stadials we find coarser material at Site U1385 than at U1391, implying that the MOW plume depth is situated between the two sites, sorting the sediment that reaches U1385 but not significantly affecting U1391 grain sizes (Figures 8c and 9). K max declinations at U1391 indicate southward flow at these times. We hypothesize that, in the transition interval, from~130 to 60 ka, the average depth of MOW deepened on the West Iberian Margin to an arrangement wherein the main plume was present at Site U1391 during interstadials ( Figure 8b) but deepened below the site during stadials (Figure 8c).

Driving Factors of the Distinct MOW Regimes
The overall structure of our records strongly suggests that the insolation response to amplitude modulation of precession by the 400 kyr eccentricity cycle exerts strong climatic control on MOW variability ( Figure 5). Longer records are needed to verify this interpretation, but MOW Regime I is closely associated with the long-eccentricity node-induced minima in modulated precession (MIS 11 and MISs 3 to 1), whereas Regime II occurs when amplitudes of the precession cycles are larger. This observation suggests a strong climate influence on MOW behavior with some kind of hydroclimate control on the salinity of the Mediterranean a likely driving mechanism. Support for this interpretation comes from the strong correspondence in structure between our records for the last glacial cycle and the reconstruction of North African rainfall climate from an eastern Mediterranean site of Ehrmann et al. (2017). In this interpretation of our data, drier conditions from about 60 ka preconditioned eastern Mediterranean waters to produce a denser deeper-equilibrating MOW in Regime I compared to Regime II. Yet, while the close correspondence between the youngest portion of our data (the transition interval and Regime I) and the records of Ehrmann et al. (2017) imply a link to North African hydroclimate and therefore summer insolation change, the strong precession signal in our records for Regime II shows a particularly close relationship with changes in spring insolation ( Figure 5). This observation raises the intriguing possibility that precipitation change associated with winter storm track activity (Bosmans et al., 2015(Bosmans et al., , 2020Kutzbach et al., 2014) may have exerted an important influence on Mediterranean hydrography and MOW behavior during the late Pleistocene.
Our data support the suggestion of Rogerson et al. (2012) that across glacial-interglacial timescales the North Atlantic vertical density structure was the primary control on the MOW plume settling depth. However, we propose that, on 400 kyr eccentricity timescales, circum-Mediterranean humidity and MOW density properties play a major role in controlling the equilibration depth of MOW in the Atlantic, via the eccentricity modulation of precession cycles.
MOW records from IODP Sites U1386 (water depth 561 m) and U1389 (water depth 644 m) in the Gulf of Cadiz (Figure 1), which extend to~570 and~250 ka, respectively, also show precession pacing and a change in behavior between~130 and 60 ka Kaboth et al., 2017). These phases and mechanisms fit well with our West Iberian Margin records of strength and depth changes of MOW. From 130 ka to present, MOW was interpreted to be more active at Site U1386, as reflected in increased grain-size variability. An offset in benthic δ 18 O values between Sites U1386 and 967 was also observed. A greater proportion of western Mediterranean deepwater (WMDW) influence, and therefore greater density contrast across the Strait of Gibraltar and faster flowing MOW, was invoked to explain these observations (Kaboth et al., 2017). During the last glacial cycle (Regime I), greater aridity in northeast Africa compared to previous glacial-interglacial cycles (Ehrmann et al., 2017;Larrasoaña et al., 2003) and the resultant higher salinity surface waters in the eastern Mediterranean would have driven more vigorous production of LIW. In this scenario, LIW would have flowed over the WMDW more quickly, causing a greater proportion of it to be drawn over the Camarinal Sill by Bernoulli aspiration (Millot, 2009;Stommel et al., 1973). This would also have acted to increase the density of the MOW, resulting in a generally lower settling depth on the west Iberian slope, as observed at Site U1391 during Regime I.
From~475 to 130 ka, Site U1386 appears to show decreased MOW activity with a greater proportion of LIW influence compared to the interval after 130 ka (Kaboth et al., 2017). This interpretation is based on both a reduced gradient between U1386 benthic δ 18 O values and those of the eastern Mediterranean ODP Site 967 and modest grain-size variability during this interval. This interpretation is broadly consistent with the shallower MOW regime that we document between~365 and 130 ka on the West Iberian Margin. Compared with Regime I, the generally more humid conditions in northeast Africa (Larrasoaña et al., 2003) during Regime II would have resulted in a less dense, and therefore slower moving, eastern sourced MOW with a much smaller proportion of WMDW, as suggested by Kaboth et al. (2017). The generally lower density MOW would have resulted in the shallower MOW settling depths that we reconstruct on the West Iberian Margin during Regime II.
Our observations call into question the interpretations of Guo et al. (2017). MOW was interpreted to be more active at Site U1391 during interglacial intervals and weaker at glacial maxima, based upon the presence of certain MOW-sensitive species of benthic foraminifera. However, a reassessment of benthic assemblages is merited (García-Gallardo et al., 2017a, 2017b because of evidence for downslope transport on the West Iberian Margin (Hall & McCave, 2000;Hodell et al., 2013;Thomson et al., 1999).
Our records clearly evidence strong climate-induced variability of MOW behavior during the late Pleistocene on both astronomical (mainly precession) and millennial timescales. Analysis of the detailed mechanisms involved is beyond the scope of this contribution, but we find strong evidence suggesting that circum-Mediterranean rainfall climate change was a likely driver of MOW variability.

Conclusions
We present new records of MOW variability on the West Iberian Margin, documenting MOW behavior back to MIS 11. The combination of AMS at Site U1391 and grain-size proxies at both Sites U1391 and U1385 allows us to reconstruct changes in MOW strength and depth. We show that AMS is a useful proxy for reconstructing current changes, mode of deposition, and migration of the MOW. Azimuth changes of the K max axis in conjunction with grain size and magnetic fabric indicate MOW presence/absence and flow speed change.
We infer two regimes of MOW behavior driven by insolation change in response to modulation of precession by the long (~400 kyr) eccentricity cycle. Strong precession-paced changes in MOW strength are documented between MISs 10 and 4 (Regime II). During this interval, interglacial peaks were marked by MOW shoaling to depths above Site U1391, indicated by the southward imbrication of magnetic grains caused by the underlying NADW at Site U1391. From MIS 5 we also see strong millennial-scale variability, which is initially superimposed on the orbital signal (defining the transition interval), but from MIS 3 to present and during MIS 11, the millennial scale signal dominates (Regime I).
Our new records signal a generally less dense more shallowly equilibrating MOW in the North Atlantic between MISs 10 and 4, compared with MISs 3 to 1 and MIS 11. Thus, the impact of MOW on North Atlantic circulation and structure may well have evolved in ways not represented during the last glacial cycle. MOW's ability to enhance NADW formation and to help reinvigorate AMOC during cold-warm transitions may have been different in Regime II. Further investigation is merited of the forcing mechanisms involved, but our records clearly document strong climate-induced variability in MOW behavior during the late Pleistocene and suggest that changes in rainfall climate in the circum-Mediterranean exerted an important control.