The impact of abrupt deglacial climate variability on productivity and upwelling on the southwestern Iberian margin

This study combines high-resolution records of nannofossil abundances, oxygen and carbon stable isotopes, core scanning X-ray fluorescence (XRF), and ice rafted debris (IRD) to assess the paleoceanographic changes that occurred during the last deglaciation on the SW Iberian Margin. Our results reveal parallel centennial-scale oscillations in coccolithophore productivity, nutricline depth and upwelling phenomena not previously observed, explained by means of arrival of iceberg-melting waters, iceberg-induced turbulent conditions, SST changes and riverine discharges. On millennial time-scales, higher primary productivity (PP), shallower nutricline, and upwelling occurrence/invigoration are observed for the Last Glacial Maximum (LGM) and Bølling-Allerød (B/A). The opposite scenario (i.e., lower productivity, deeper nutricline and upwelling weakening/absence) is linked to cold spells such as Heinrich Stadials 2 and 1 (HS2 and HS1) and the Younger Dryas (YD). Such paleoproductivity variations are attributed to latitudinal migrations of the thermal fronts associated with oceanic gyres in the North Atlantic, in parallel to oscillations in the strength of the Atlantic Meridional Overturning Circulation (AMOC). Moderate-to-high PP during the Holocene is ascribed to the development of the modern seasonal surface hydrography, with a more persistent Iberian Poleward Current (IPC) and seasonal wind-induced upwelling. © 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Marine sedimentary sequences recovered from the SW Iberian margin have played a pivotal role in unraveling the pace, duration and magnitude of past climate changes (e.g., Abrantes et al., 2010;Hodell et al., 2013;Martrat et al., 2007;Shackleton et al., 2000). These sediments have been explored extensively to decipher past environmental conditions that occurred during abrupt (millennial) climate events from the LGM thorough the last deglaciation (e.g., Abrantes et al., 2001;Boessenkool et al., 2001;Eynaud et al., 2009;Naughton et al., 2016;Salgueiro et al., 2014;S anchez-Goñi et al., 1999;Sch€ onfeld and Zahn, 2000). However, little is known about the impact of abrupt climate events and related environmental conditions on the local upwelling system, for which highresolution productivity records are lacking and conflicting interpretations exist for specific time periods. For instance, during Heinrich Stadials (HSs), high primary productivity (PP) and export productivity (Pexp) have been interpreted from proxy records, ascribed to turbulent mixing resulting from iceberg drift, the southward penetration of a nutrient-rich water mass and/or nutrient input from the shelf or fluvial discharge (Eynaud et al., 2000;Lebreiro et al., 1997;Thomson et al., 2000). Conversely, for the same period, other proxy records have been interpreted as a reduction in PP and Pexp, as a consequence of upwelling cessation linked to strong surface water stratification (Incarbona et al., 2010;Pailler and Bard, 2002;Salgueiro et al., 2010;Voelker et al., 2009). Additionally, Salgueiro et al. (2014) suggested the presence of frontal upwelling to explain latitudinal differences in Pexp along the Iberian margin during HSs. The Last Glacial Maximum (LGM), a period of more moderate stadial conditions compared to HS2 preceding and HS1 proceeding it, is characterized by higher PP attributed to intensified trade winds and subsequent upwelling (Abrantes, 1991;Palumbo et al., 2013;Voelker et al., 2009). Few PP and Pexp reconstructions are available for the Younger Dryas (YD), although existing Iberian Margin records demonstrate contrasting estimations (Incarbona et al., 2010;Palumbo et al., 2013;Salgueiro et al., 2014), which are typically not sufficiently resolved to confidently address productivity variations during its several phases (Naughton et al., 2019). Despite the majority of paleoproductivity reconstructions for the last deglaciation focus on the stadials, high PP is observed for the interstadial Bølling/Allerød (B/A) (Incarbona et al., 2010;Pailler and Bard, 2002;Palumbo et al., 2013). Incarbona et al. (2010) hypothesized two non-excluding scenarios to explain increased PP during glacial interstadials: 1) eastward migration of the Azores Current (AzC) promoting frontal upwelling and enhancing PP and 2) reinforcement of the northern Hemisphere atmospheric circulation leading to a widening of the upwelling cell near Cape St. Vincent (Fig. 1). Finally, the Holocene features decreased (Palumbo et al., 2013), moderate (Salgueiro et al., 2010;Voelker et al., 2009), and variable upwelling and productivity levels, the latter predominately consisting of oligotrophy during the Early Holocene, with increased PP during the Late Holocene (Palumbo et al., 2013(Palumbo et al., , 2019Rodrigues et al., 2009).
Each of these abrupt climate changes has been associated with shifts in the strength of the Atlantic Meridional Overturning Circulation (AMOC) in the North Atlantic (McManus et al., 2004;Ng et al., 2018), including the Iberian margin (Gherardi et al., 2005). Proxy evidence and model simulations suggest the weakening and/ or shutdown of AMOC via freshwater input in the North Atlantic during HS2, HS1 and the YD led to Northern Hemispheric cooling, decreased precipitation and reorganization of the trade wind circulation (Jackson et al., 2015;Menviel et al., 2008;Naughton et al., 2016Naughton et al., , 2019. In the ocean, AMOC reduction has been linked to latitudinal migrations of North Atlantic oceanic gyres (Reibig et al., 2019;Repschl€ ager et al., 2015) and a large reduction in export production, the latter ascribed to shoaling of the mixed layers and the subsequent decreased in nutrient input (Menviel et al., 2008;Schmittner, 2005).
As one of the main phytoplanktonic groups, coccolithophores can be used to decipher upwelling strength and primary productivity variations . On the modern-day W Iberian margin, these calcifying algae bloom during the final stage of the upwelling event, when warming and stratification of surface waters begin and coccolithophores outcompete diatoms for the remaining nutrients Silva et al., 2008Silva et al., , 2009. In contrast, nutrient scarcity and strong vertical mixing during the downwelling regime hamper coccolithophore proliferation. Thus, higher accumulation rates of calcareous coccolithophore remains (referred to as calcareous nannofossils) can be interpreted in the sedimentary record as periods of increased PP (Amore et al., 2012;Palumbo et al., 2013). Previous investigations using sediment trap, water column, and surface sediment samples have assessed the ecology of extant coccolithophore species in the study area Cachão and Moita, 2000;Ferreira and Cachão, 2005;Guerreiro et al., 2013Guerreiro et al., , 2014Moita et al., 2010;Narciso et al., 2006;Silva et al., 2008). Information gleaned from the latter studies facilitate the interpretation of variations in the abundance of key species in the fossil record in terms of the specific oceanographic and climatic conditions that preferentially favor or hamper their proliferation.
We examined sediments from core SHAK06e5K from the SW Iberian margin, whose chronostratigraphy provides accurate and robust age control on the derived proxy records (Ausín et al., 2019a). By combining high-resolution nannofossil records with other sedimentological and inorganic geochemical data (carbon and oxygen isotopes, XRF and IRD counts) from the same sediment core, we increased the resolution of previous paleoproductivity reconstructions and document the impact of rapid climate transitions on PP and surface water dynamics on the SW Iberian margin over the last 28 kyr.

Hydrographic conditions
The SW Iberian margin is located on the north-eastern edge of the North Atlantic subtropical gyre (Fig. 1). Surface waters in the region are bathed by the Eastern North Atlantic Central Water (ENACW) (van Aken, 2001). Throughout the year, the Portuguese Current (PC), which branches from the relatively cold North Atlantic Current, flows southward (P erez et al., 2001). From October to March, the Iberian Poleward Current (IPC), which is fed by the eastern branch of the warmer AzC, flows northward close to the coast along the W Iberian margin (Haynes and Barton, 1990). The region is influenced by two recurrent thermal frontal systems (Peliz et al., 2005): the Subtropical Front, which is associated to the eastern branch of the AzC and fluctuates along 35e36 N, and the Western Iberian Winter Front, which separates cold waters along 40-38 N. During times of IPC circulation, the W Iberian Winter Front deflects northward along the Iberian margin separating this current from the PC (Peliz et al., 2005).
As part of the Canary Eastern Boundary Upwelling Ecosystem, one of the world ocean's four major coastal upwelling systems, the study area sustains high productivity (Arístegui et al., 2009). Current PP is driven by seasonal changes in the atmospheric system. From March/April to September/October, wind-induced upwelling of sub-surface waters dominate the hydrography, fertilizing the surface ocean and enhancing PP (Relvas et al., 2007). From September/October to March/April, prevailing southwesterly winds promote coastal downwelling (Peliz et al., 2005).

Age-depth model
The age-depth model for SHAK06e5K is a depositional model (P_sequence) based on 41 AMS 14 C dates of monospecific samples of Globigerina bulloides and built with the calibration package Oxcal (Bronk Ramsey, 2009). Conventional radiocarbon ages were calibrated using the Marine13 curve, which applies a marine reservoir correction of 400 yr (Reimer et al., 2013), while no local reservoir was applied. Conventional radiocarbon and calibrated ages have been published elsewhere (Ausín et al., 2019a(Ausín et al., , 2019b. According to Skinner et al. (2014) marine reservoir ages on the Iberian margin might have been greater for HS1 and YD, but here, similar to the vast majority of previous works focusing on the study region, a static correction has been applied to favor comparison with nearby records. Accordingly, discrepancies are expected between the calendar ages estimated for HS1 and YD in our records and other published works where different assumptions about the marine reservoir correction have been made (i.e., Hodell et al., 2017). The age-depth model spans the last 28 kyr and resulting sedimentation rates vary between 10 and 30 cm/kyr for most of the core, decreasing to 3e6 cm/kyr during the middle and late Holocene. A 1-cm sampling resolution yields an average temporal resolution of 120 yr for the bottom core (from 328 to 88 cm) and 310 yr for the most recent section (from 88 to 0 cm).

Nannofossil-based productivity and environmental proxies
A total of 156 samples were studied for calcareous nannofossil analyses. Samples were prepared following the settling technique of Flores and Sierro, 1997. Counts were performed with a Zeiss Axio polarizing microscope with a phase contrast device at 1000Â magnification. To ensure a confidence level of 99% for species with a proportion of at least 2% of the assemblage, a minimum of 450 specimens per sample were counted and identified, then used to calculate relative (%) and absolute (coccoliths g À1 ) abundance of each taxa. Relative abundance is exclusively interpreted as an indicator of assemblage composition at a given time. Nannofossil accumulation rate (NAR, given in coccoliths cm À2 ka À1 ) was calculated by considering the absolute abundance of each taxa, sedimentation rate and dry bulk density of the sediment.
Species identification follows the taxonomic criteria of Young et al. (2003). Following other studies (Amore et al., 2012;Palumbo et al., 2013), total NAR is used as indicator of coccolithophore productivity and abundance variability through time. The study of monthly time-series of water samples from the Iberian margin unequivocally linked Emiliania huxleyi to the upwelling regime (Silva et al., 2008), particularly at the end of upwelling events, when the water column starts to warm and stratify . Accordingly, E. huxleyi is used as an upwelling indicator. G. aperta and G. ericsonii were combined within the small Gephyrocapsa group (<3 mm). These taxa were linked to eutrophic conditions related to enhanced upwelling , and are typically used as an indicator of nutrient-rich waters in the region (Amore et al., 2012;Colmenero-Hidalgo et al., 2004;Palumbo et al., 2013). Similarly, Coccolithus pelagicus ssp. braarudii is used as a proxy of upwelled waters in the study area (Amore et al., 2012;Cachão and Moita, 2000;Parente et al., 2004;Silva et al., 2008). Ausín et al. (2018) related G. oceanica to nutrient content and water depths below 50 m. This species has typically been used as a paleoproductivity proxy (Amore et al., 2012, and references therein). Florisphaera profunda is a lower photic zone dweller whose abundance in water samples from the NW Iberian margin has been unequivocally linked to the downwelling regime and low-productivity conditions . A recent global calibration model, based on core-top samples, allows quantitative estimation of net PP from F. profunda percentages in fossil cores (Hern andez-Almeida et al., 2019). However, the optimal inverse relationship between F. profunda and net PP is limited to the latitudinal range between 30 S and 30 N, and thus, its relative abundance in our records is only interpreted as a qualitative indication of nutricline depth and PP (Beaufort and Heussner, 2001;Incarbona et al., 2010, and references therein;Marino et al., 2014). Gephyrocapsa muellerare is a nutrient-rich, cold-water indicator (Palumbo et al., 2013;Silva et al., 2008;Weaver and Pujol, 1988), although its abundance varies regardless of absolute SST changes (Incarbona et al., 2010). E. huxleyi large (>4 mm) has been linked to both, cold (Colmenero-Hidalgo et al., 2002 and cold and fresher (Ausín et al., 2015;Bazzicalupo et al., 2018) waters, whereas Coccolithus pelagicus ssp. pelagicus is a typical cold water proxy (Amore et al., 2012;Colmenero-Hidalgo et al., 2004;Palumbo et al., 2013;Parente et al., 2004). The taxa Umbilicosphaera spp., Umbellosphaera spp., Syracosphaera pulchra, Oolithotus spp. and Rhabdosphaera spp., were grouped within the "warm water group" (WWG) according to their common affinity for waters of subtropical origin .
The N ratio was calculated according to Flores et al. (2000): Firstly, it considers the abundance of the small placoliths (i.e. E. huxleyi and small Gephyrocapsa), which are typical upper photic zone (UPZ) dwellers characteristic of high-productivity conditions. Secondly, it reflects the abundance of the lower photic zone dweller F. profunda that occupies a deep nutri-thermocline position. Thus, higher N ratio values imply a deeper nutricline.
Specimens obviously pertaining to older stratigraphic levels are referred to as "reworked specimens". Their relative abundance in relation to the autochthonous assemblage is used as an indicator of particle advection to the study site (Ferreira et al., 2008;Incarbona et al., 2010).

XRF core scanning
Core SHAK06-5K was subsampled along its length (344 cm) with three u-channels (plastic core-length boxes with 2 cm Â 2 cm cross sections and maximum 150 cm length) onboard ship. Each uchannel surface was carefully scraped cleaned and covered with a 4 mm thin SPEXCertiPrep Ultralene foil to avoid contamination and minimize desiccation. Semi-quantitative elemental data were obtained using an Avaatech X-ray fluorescence (XRF) core scanner at the Godwin Laboratory in Cambridge (U.K.). Each section was irradiated with a rhodium X-ray source using a tube current of 0.2 mA at three different voltages: 10 kV (kV) with no filter, 30 kV using a thin lead filter, and 50 kV using a copper filter. XRF data were collected every 5 mm along the entire length of each u-channel. The length and width of the irradiated surface was 2.5 mm (down-core) and 12 mm (cross-core), respectively, with a count time of 60 s. The Canberra WinAxil software was used with standard software settings and spectrum-fit models to obtain element intensities, which depend on the element concentrations. Results are presented in log ratios of element intensities, which most accurately reflects changes in chemical composition (Weltje and Tjallingii, 2008).

IRD counts
Ice Rafted Debris (IRD) counts were performed using a light microscope at 2-cm resolution in the sediment fraction >300 mm from the depth intervals encompassing HS2 (308-262 cm), HS1 (170-114 cm), and YD (96-78 cm). Detrital carbonate grains were counted as a sub-category of the total lithic grains. None of the samples were split. IRD concentration is expressed in grains per dry fraction weight (grains g À1 ; >300 mm).

Oxygen and carbon stable isotope records
High-resolution carbon (d 13 C) and oxygen (d 18 O) stable isotope records of the planktonic foraminifer G. bulloides from core SHAK06e5K have been presented elsewhere along with details on sample preparation and analyses (Ausín et al., 2019a). Both records reveal centennial-to millennial-scale oscillations linked to the major abrupt climate events of the last deglaciation (eg., HS2 and HS1, B/A, and YD, respectively). The planktonic d 18 O variations reflect changes in sea surface temperature (SST), ice-volume and salinity. Planktonic foraminifera d 13 C is assumed to reflect the carbon isotopic composition of dissolved inorganic carbon (d 13 C DIC ) as well as 'vital' effects, which include kinetic isotope fractionation that results in isotopic disequilibrium (Spero, 1992). Because d 13 C DIC depends on local (i.e. photosynthesis and respiration, changes in advection and upwelling, calcification rate) and global carbon cycle changes, d 13 C shell reflects multiple parameters (Ravelo and Hillaire-Marcel, 2007). At given locations where upwelling strongly impacts surface local hydrography, d 13 C is often interpreted as a proxy of upwelling (Prasanna et al., 2016;Voelker et al., 2009;Zhang et al., 2017). Accordingly, we use d 13 C, along with other paleoproductivity and nutricline depth indicators, to interpret variations in the upwelling dynamics.

Nannofossil abundances
Preservation of specimens was good to moderate according to the visual criteria established by Flores and Marino (2002). Relative abundance and NAR of specific taxa used as paleoceanographic proxies in this study (see section 3.2.1) are shown in Figs. 2 and 3. E. huxleyi and the small Gephyrocapsa group are the most abundant taxa throughout the studied period, comprising up to 80% and 50% of the assemblage, respectively (Fig. 2).
The LGM is again dominated by the small placoliths of E. huxleyi and small Gephyrocapsa (50% and 30% on average, respectively), but  the assemblage is variable throughout this period. For instance, larger percentage contribution of Coccolithus pelagicus ssp. braarudii (~1.4%) and G. oceanica (~5%) are observed from 23.5 to 21.5 ka (Fig. 2). Between 21.5 and 19.5 ka, only Coccolithus pelagicus ssp. braarudii exhibits a clear increase in relative abundance (up to 2%).
From 19.5 to 18.5 kyr, G. muellerare, E. huxleyi (>4 mm) and G. oceanica show higher NAR and percentages, whereas smalls peaks of WWG are observed throughout the entire LGM. Peaks up to 1.8% of Coccolithus pelagicus ssp. braarudii are observed during the B/A, a period marked by a higher relative abundance of small Gephyrocapsa (up to 35%) and the WWG (~4%). The first half of the YD (YDa; 13.2 to 12.5 ka) is characterized by a decrease in the relative abundance of most species and a large increase in E. huxleyi, which reaches up to 70% at 12.7 ka (Fig. 2). G. oceanica and F. profunda feature an increasing trend during the second half of the YD (YDb; 12.5 to 11.4 ka), comprising of 4% and 10%, respectively, of the assemblage at the end of the interval. The WWG constitutes 4% of the assemblage during YDb (Fig. 3). The Holocene assemblage is mainly composed of E. huxleyi (70% on average) with lesser but significant contributions of F. profunda (20%), WWG (5% on average), and G. oceanica (3%). The relative abundance of reworked specimens has been presented elsewhere (Magill et al., 2018). Higher percentages are observed from 28 to 12.5 ka, especially during HS2 (~10%) and HS1 (~16%) and the first half of the YD (4%), and decrease to minimum from 12.5 ka onwards (Fig. 3).
In general, total NAR displays higher values during the LGM, B/A and most of the late Holocene (last 3.7 kyr) (Fig. 4c). Three larger peaks are observed at 22.1 ka, 20.7 ka, and 19 ka for the LGM. The Holocene is marked by low values until 3.7 ka, from which high variability and larger peaks are observed until present. NAR of the majority of the taxa follow this same pattern, with the following exceptions: G. muellerare and C. pelagicus ssp. pelagicus are generally more abundant from 28 to 12.5 ka and show little increase from 3.7 ka onward, while F. profunda and the WWG also show a significant increase in NAR along the Holocene in relation to the deglaciation (Figs. 2 and 3). This time span is characterized by the near absence of G. muellerae and E. huxleyi (>4 mm) (Fig. 3), and lower NAR of small Gephyrocapsa group, C. pelagicus ssp. braarudii (Fig. 2) and C. pelagicus spp. pelagicus (Fig. 3). Other taxa such as E. huxleyi, G. oceanica, F. profunda and the WWG increase their NAR and variability during the last 3.7 kyr.
The N ratio exhibits a double-peak decrease during both HSs, which is more pronounced during HS1 (Fig. 4d). Maximum values are observed during the LGM, B/A and YDa. A general decrease in N ratio is observed from the onset of YDb thorough the Holocene, reaching a minimum at 3 ka with a relative recovery to high values thereafter.

XRF and IRD counts
The log (Ca/Ti) was selected to identify variations in %CaCO 3 (Hodell et al., 2013(Hodell et al., , 2015, whereas log (Ca/Sr) allows the relative proportion of detrital material (dolomite and limestone) to be distinguished from biogenic carbonate (Hodell and Curtis, 2008;Hodell et al., 2017). The Ca/Ti record shows a general increasing trend across the last deglaciation with the lowest values occurring in early HS1, followed by a peak in the middle of HS1 (~17 ka), decreasing again during the latter part of HS1 (Fig. 4g). Ca/Ti values increase during the B/A, followed by a decline during the YD before peaking in the early Holocene. In contrast, Ca/Sr is relatively low throughout the record with the exception of a prominent increase centered at 17 ka (Fig. 4h). Both Ca/Sr and Ca/Ti records show a conspicuous peak during HS1 at 17 ka. Ca/Sr also shows a slight peak during HS2 at 24.5 ka that is not expressed in the Ca/Ti signal.
Lithic grain counts were targeted in the intervals suggestive of IRD delivery on the basis of XRF results. These results demonstrate significant amounts of IRD from 25.2 to 24.3 ka during HS2 and 17.5 to 16 ka during HS1, reaching absolute maximum values at 17 ka and 25 ka, respectively (Fig. 4i). No IRD was found during the YD event.
5. Discussion 5.1. Heinrich stadials 2 and 1 5.1.1. General conditions HS2 and HS1 are characterized by minimum coccolithophore productivity, a deeper nutricline (Fig. 4), and a general decrease in Ca/Ti (Fig. 4) indicating a larger input of detrital sediment (Hodell et al., 2013(Hodell et al., , 2015. Excluding additional effects, significant increases in d 13 C of G. bulloides are interpreted as nutrient scarcity in the upper water column and hence decreased upwelling of deeper waters (Ausín et al., 2019a;Voelker et al., 2009). Such conditions were likely to hinder not only coccolithophore productivity, but overall PP and Pexp, as inferred by previous studies in this region (Incarbona et al., 2010;Salgueiro et al., 2010;Voelker et al., 2009) and other northern locations of the North Atlantic (Nave et al., 2007). It is likely that variations in the upwelling system were primarily induced by oceanic (i.e, hydrologic changes) rather than atmospheric forcing (e.g., pressure center dynamics and subsequent blowing winds). Hydrologic changes in the physical structure of the upper water column, such as strong stratification, have been inferred from the arrival of cold, less-saline waters . Melting icebergs derived from the northern ice-sheets delivered meltwater and ice rafted debris to the study region (Bard et al., 2000;Cayre et al., 1999;Eynaud et al., 2009;Voelker and de Abreu, 2011). Such conditions are corroborated by the presence of IRD coinciding with d 18 O decreases for both HS2 and HS1 (Fig. 4). According to Ferreira et al. (2008) and previous work undertaken in core SHAK06e5K (Magill et al., 2018), the larger contributions of reworked specimens observed at these times are the result of increased lateral transport of fine-silt size (<10 mm) sediments in relation to intermediate nepheloid layers associated with increased MOW transport during HS2 and HS1 (Fig. 3).

Differing impact on productivity and hydrography between HS2 and HS1
During HS1, a deeper nutricline position and reduced/suppressed upwelling is demonstrated in comparison to HS2 (Fig. 4). These conditions result from the larger impact of HS1 on the regional hydrography, as indicated by greater IRD abundances and a larger d 18 O decrease due to a more significant freshwater input, which point to a greater intensity and duration of HS1 compared to HS2. Similarly, significant delivery of detrital carbonate d a diagnostic feature of Heinrich layers (Hodell and Curtis, 2008;Hodell et al., 2017) d is only visible during HS1. Such conditions induced substantial differences in nannofossil assemblages between both HS events. For instance, small Gephyrocapsa and G. muellerae largely contribute to the assemblage during HS2, indicating more eutrophic conditions towards the end of the period. In contrast, high percentages of E. huxleyi (>4 mm) and C. pelagicus ssp. pelagicus are only observed during HS1. Both species are indicators of subpolar water and are known to peak during Heinrich events off the Iberian margin (Narciso et al., 2006;Parente et al., 2004) supporting a greater effect of iceberg melting at this time accompanied by a larger SST decrease (Salgueiro et al., 2014). These features, along with a major southward migration of the polar front during HS1 (Eynaud et al., 2009) may have posed an ecological barrier for the development of certain species.

Centennial productivity variations during HSs
Some of the proxy-records exhibit a bimodal distribution during both HSs that deserves further attention. Double peaks of F. profunda percentages led to decreases in the N ratio, coinciding with higher d 13 C excursions centered at 25.6 and 24.7 ka for HS2, and at 18 and 16 ka for HS1 (Figs. 2 and 4d and e), indicating large centennial variations in upwelling intensity and nutricline depth. The intervening intervals are characterized by the arrival of iceberg and meltwater, indicated by IRD and XRF peaks and a slight decrease in d 18 O (Fig. 4fei). The IRD and XRF peaks at 17 ka ( 14 C age of 14,089 ± 101) likely correspond to those identified in the North Atlantic at 16 ka cal. BP (corresponding 14 C age of 13,913 ± 30) by Hodell et al. (2017, and references therein) and termed H1.1. The age discrepancy in the calendar ages result from the different reservoir ages applied in each case (1200 yr in Hodell et al. (2017) and a static reservoir of 400 yr in this study). Double peaks in upwelling weakening and nutricline deepening are attributed to the influence of meltwaters right before and after maximum iceberg presence. Fresher surface waters are expected to promote a stronger halocline, suppress wind-induced upwelling, and inhibit PP and Pexp (Incarbona et al., 2010;Pailler and Bard, 2002;Salgueiro et al., 2010;Voelker et al., 2009), which is consistent with the overall low coccolithophore productivity observed in our records (Fig. 4c). During maximum iceberg intrusion (25 and 17 ka), however, a relative nutricline shoaling and upwelling intensification is observed (dark blue bars in Fig. 4d and e), in contrast with the lowest coccolithophore productivity observed in our records (Fig. 4c). Turbulent mixing due to iceberg drift (Eynaud et al., 2000;Lebreiro et al., 1997;Thomson et al., 2000) might account for such variations in the nutricline depth and upwelling intensity, whereas it is likely that overall coccolithophore productivity was hampered by the arrival of colder and fresher waters, in which only E. huxleyi (>4 mm) and C. pelagicus ssp. pelagicus would find optimal growth conditions ( Fig. 3). Given that certain primary producers (e.g. dinoflagellates and diatoms) can outcompete coccolithophores under medium and high turbulent conditions (Villamaña et al., 2019), the proposed scenario allows the reconciling of previous contrasting productivity interpretations for the study area: low PP and Pexp inferred from coccolithophore, planktonic foraminifera, alkenone accumulation and CaCO 3 data (Incarbona et al., 2010;Pailler and Bard, 2002;Salgueiro et al., 2010;Voelker et al., 2009), and higher PP and Pexp interpreted from dinoflagellate cysts and planktonic foraminifera species (Eynaud et al., 2000;Lebreiro et al., 1997;Salgueiro et al., 2014). Parallel centennial oscillations in reworked specimens and WWG also deserve further attention. Deeper position of the MOW flow core, deduced from percentage increases in reworked specimens in core SHAK06-5K (Magill et al., 2018), may have enabled the eastward displacement of AzC (Incarbona et al., 2010, and references therein). Discrete peaks in WWG at these times would then be explained by the short-lasting influence of oligotrophic and warmer subtropical waters.

The LGM
Coccolithophore productivity increased significantly during the LGM, linked to a shallow nutricline and an overall intensification of upwelling (Fig. 4cee). Such conditions suggest an increase in marine productivity off W Iberia, as previously shown from both coccolithophore (Palumbo et al., 2013) and diatom fossil records (Abrantes, 2000), as well as from Pexp estimated from a planktonic foraminifera-based transfer function (Salgueiro et al., 2010;Voelker et al., 2009). Relatively stable and mild SSTs prevailed during this period (Darfeuil et al., 2016), albeit with temperatures low enough to promote proliferation of cold-water species. Coccolithophore productivity and upwelling intensity show large centennialmillennial oscillations that are superimposed on these general trends. For instance, NAR shows three distinct peaks centered at 22.1 ka, 20.7 ka, and 19 ka, which are the result of large changes in sedimentation rate along with concomitant increases in coccolith abundance (coccolith g À1 ) (Fig. 4aec). Each of these peaks is seen in the NAR of E. huxleyi and small Gephyrocapsa, suggesting high productivity at those times. Nevertheless, changes in the assemblage also indicate temporally varying environmental conditions. Small peaks of WWG along the entire period indicate the seasonal incursion of warmer, subtropical waters, possibly related to a Paleo-IPC (Salgueiro et al., 2014;Voelker et al., 2009). More eutrophic conditions, inferred from the larger percentage contributions of small Gephyrocapsa, G. oceanica, and Coccolithus pelagicus ssp. braarudii, are maintained for the first part of the period until 21 ka (Fig. 2). At this time, a transition to colder waters is inferred from the large increase in NAR and relative abundance of classical coldwater indicators G. muellerare and C. pelagicus ssp. pelagicus (Fig. 3).
The major increase in E. huxleyi (>4 mm) abundance commencing at 19 ka indicates the influence of not only cold but fresher conditions (Fig. 3). Concomitant peaks of this taxa in Western Mediterranean cores (Ausín et al., 2015;Bazzicalupo et al., 2018;Flores et al., 2010a) have been linked to major fluvial discharge from North Europe glaciers (Bazzicalupo et al., 2018). Our high-resolution record supports massive glacier runoff from the British Islands and Scandinavia into the North Atlantic (Bazzicalupo et al., 2018), flowing down to the study area and entering the Mediterranean through the Strait of Gibraltar.

The Bølling/Allerød
The B/A was marked by high coccolithophore productivity and a shallow nutricline linked to increased upwelling (Fig. 4cee). More eutrophic conditions in relation to anterior HS1 are supported by increases in the relative abundance of small Gephyrocapsa and C. pelagicus ssp. braarudii, in agreement with previous PP qualitative estimations (Incarbona et al., 2010). Continuous SST increase from the Bølling to the Allerød has been observed from sediment cores along the W Iberian margin (e.g., Martrat et al., 2014;Naughton et al., 2016;Rodrigues et al., 2010). This substantial warming is accompanied by a gradual transition from subpolar (C. pelagicus ssp. pelagicus) to subtropical (WWG) taxa (Fig. 3), likely promoted by the northward extension of the subtropical gyre (Schwab et al., 2012). The subsequent northward displacement of the associated Azores Front would explain the continuous increase in PP and nutricline shallowing observed along the B/A.

The YD
The YD is a major cooling event in the northern Hemisphere (Alley, 2000) marked by a 4 C drop in SST in the study area (Darfeuil et al., 2016). Overall productivity is low, accompanied by ongoing nutricline deepening and reduced upwelling ( Fig. 4c and  d). Conversely, Palumbo et al. (2013) observed a conspicuous increase in NAR for the entire YD from sediment core MD03-2699 ( Fig. 1), inferring increased PP. Similar latitudinal differences have been recognized in the quantitative reconstruction of Pexp in northern and southern sites off the Portuguese coast Salgueiro et al., 2014). It has been suggested that the polar front was not displaced as far south as our study site (Eynaud et al., 2009;Naughton et al., 2016). Yet, the arrival of colder and less saline waters (Duplessy et al., 1992) affected the W Iberian margin from north to south. These conditions likely promoted increased surface stratification rather than upwelling of nutrient-rich waters. Nevertheless, a recent study based on pollen records reveal asymmetrical north/south rainfall patterns along the margin (Naughton et al., 2016). It is possible that the extreme winter precipitation that affected the north and central parts of the margin led to strong coastal nutrient input by fluvial discharges (Rodrigues et al., 2010), leading to increased productivity in those regions despite surface stratification and reduced/absent upwelling, contrasting with the general low productivity in SW Iberia.
In greater detail, our records enable at least two distinct phases during the YD to be distinguished, named here YDa (13.4e12.5 ka) and YDb (12.5e11.4 ka). YDa is the coldest phase as indicated by heavier d 18 O, peaks in NAR and percentages of G. muellerae, a drop in WWG abundances, and, similar to HS2 and HS1, the relative increase in reworked coccoliths which suggests a stronger MOW. The opposite is true for YDb, suggesting warmer conditions during the second phase. Recently, the YD has been described as a four-phase interval from alkenone-derived SST and pollen records at decadal resolution from nearby core D13882 (Naughton et al., 2019). Briefly, a cool and dry onset (12890-12720 yr BP) followed by an exacerbation of these conditions (12720-12390 yr BP) constitute the first two phases, which likely correspond to YDa in this work. The third phase in Naughton et al. (2019), consisting of gradual warming and moisture increase (2390e12030 yr BP) corresponds to YDb, whereas a complex fourth phase featuring rapid changes between cool-wet, cold-dry and cool-wet conditions (12030-11770 yr BP) is not observed in our records due to its lower resolution. Despite weaker upwelling during YDa in relation to YDb, the former shows higher productivity and a shallower nutricline, in contrast with general stratified conditions. Assemblage composition suggests higher total NAR during YDa might be due to occasional blooms of the opportunistic E. huxleyi. However, higher-resolution records are needed to provide reliable productivity variations across the different phases of the YD.

The Holocene
Both the Early and mid-Holocene are marked by moderate/low coccolithophore productivity and a deeper position of the nutricline (Fig. 4). Both time periods include the warmest SST for the last 28 kyr (Darfeuil et al., 2016;Rodrigues et al., 2009). This, along with weakening of upwelling-favorable winds over Iberia (L ezine and Den efle, 1997), may have supported considerable stratification, limiting the extension of the Ekman layer and preventing strong upwelling episodes. The 2 C cooling initiated during the mid-Holocene (Darfeuil et al., 2016) possibly favored a less stratified water column and increased PP during the Late Holocene in relation to the Early and mid-Holocene (Fig. 4cee) (Palumbo et al., 2019;Rodrigues et al., 2009). However, the relatively low total NAR during the Early and mid-Holocene stems largely from the significant reduction in sedimentation rate during this time period (Fig. 4a) masking the contribution of coccolithophore absolute abundances (Fig. 4b). The latter record indicates relatively higher absolute abundances during the Early Holocene and the second half of the mid-Holocene, varying in parallel to nutricline shallowing and upwelling intensity. Superimposed onto these general conditions are striking short-term in-phase variations between these three parameters, with synchronous peaks shown at 9.9, 8.8, 3.7 and 2.5 ka, implying centennial co-variation between absolute coccolithophore abundances, nutricline depth and upwelling strength for the entire Holocene. Increased WWG abundances during the Holocene indicate the IPC became a more common feature of the seasonal surface circulation, in agreement with the northward recirculation of the AzC proposed for this period (Palumbo et al., 2013(Palumbo et al., , 2019, progressively leading to the modern hydrographic configuration, where seasonal wind-induced upwelling is the major control of PP variations in the study area.

Drivers for paleoproductivity variations on millennial-time scale
Despite the influence of different climate forcings during each of the abrupt climate events of the last glacial, a pattern arises regarding productivity variations: the lowest PP levels were found during the coldest glacial stadials (i.e., HS2, HS1 and YD), whereas the warmer LGM and interstadial B/A were characterized by higher PP. This pattern is in line with results from cores MD01-2444 (Incarbona et al., 2010), MD95-2040and MD9-2042(Pailler and Bard, 2002 (Fig. 1), and are further supported in our records by in-phase nutricline deepening (shallowing) and upwelling absence/weakening (occurrence/invigoration) during HSs and the YD (LGM and interstadial B/A). Our results suggest that surface water-column stratification induced by iceberg melting waters account for the lowest PP observed during the coldest glacial stadials in the SW Iberian margin (Incarbona et al., 2010;Pailler and Bard, 2002;Salgueiro et al., 2010;Voelker et al., 2009), as observed for northern Atlantic sites (Nave et al., 2007). The reversed productivity pattern is observed in open sea sites under the current influence of the North Atlantic Subtropical gyre, where the northward expansion of the latter and the influence of oligotrophic waters promoted decreased PP during the B/A (Schwab et al., 2012) and interglacials (Nave et al., 2019). The subsequent northward and eastward displacement of the associated Subtropical Front would have possibly reached our study site, enabling frontal upwelling and enhanced PP during warmer periods (i.e., LGM and B/A), in line with the eastward migration of the AzC hypothesized by Incarbona et al. (2010) for interglacials. This "seesaw" in productivity patterns on millennial time-scales in the North Atlantic suggests the influence of a common forcing mechanism operating throughout the last glacial. Latitudinal migrations of the thermal fronts associated with the expansion/contraction of oceanic gyres have been ascribed to oscillations in the strength of AMOC triggered by North Atlantic iceberg discharges (Reibig et al., 2019;Repschl€ ager et al., 2015). During weak AMOC phases (i.e., HSs and YD) (Gherardi et al., 2005;Henry et al., 2016;McManus et al., 2004;Ritz et al., 2013), both the Subpolar and Subtropical fronts migrated southward and iceberg melting waters reached the Iberian margin, whereas a northward migration of both fronts occurred during strong AMOC phases (i.e., LGM and B/A) (Eynaud et al., 2009;Reibig et al., 2019;Repschl€ ager et al., 2015). Marine productivity decreases during weak AMOC phases estimated from model simulations have been attributed to the shoaling of mixed layers, the formation of a strong halocline and the subsequent reduced nutrient supply to the photic zone (Menviel et al., 2008;Schmittner, 2005). The strong coupling between latitudinal displacements of thermal fronts and opposite PP patterns in the study area and the subtropical North Atlantic (Nave et al., 2019;Schwab et al., 2012) supports a non-excluding mechanism by which AMOC would modulate PP on millennial time-scales in the North Atlantic via expansion and contraction of North Atlantic oceanic gyres.

Conclusion
Overall, a remarkably good agreement on centennial and millennial time-scales is observed between higher (lower) coccolithophore productivity, shallower (deeper) nutricline, and stronger (weaker) upwelling events for the last 28 kyr. The high-resolution records presented here allows centennial PP changes not previously observed in the study area to be identified, and demonstrates that these variations were primarily controlled by the hydrographic and environmental conditions resulting from abrupt deglacial climate events. Such short-term variations are explained by means of intrusions of iceberg-melting waters, turbulence promoted by iceberg drift, changes in SST and fluvial discharges. An obvious pattern emerges regarding PP variations on millennial time-scales: higher PP during warmer LGM and B/A and lower PP during the coldest glacial stadials (HS2, HS1 and YD). We attribute this pattern to changes in sea-ice extent and latitudinal migrations of the thermal fronts associated to the expansion and contraction of oceanic gyres, which points to a possible linkage between AMOC variability and PP changes on millennial time-scales in the study area.