North Iberian temperature and rainfall seasonality over the Younger Dryas and Holocene

Several stalagmite records have yielded important but discontinuous insights into northern Iberian climate since the Last Glacial. Here we present the first continuous Iberian stalagmite-based reconstruction of climate since the Bølling-Allerød interstadial, from a single stalagmite sample (GAR-01 from La Garma Cave, Cantabria). The ~13.5 ka GAR-01 record provides the opportunity for replication, continuation, and aggregation of previously published records from northern Spain. The GAR-01 record reveals shifts in oxygen isotope ratios that are inexplicable by appealing to a single control (i.e., exclusively temperature, rainfall amount, etc.). Herein we explore the potential role of rainfall and temperature seasonality shifts on the new d18O record using a simple Monte Carlo approach to estimate the seasonal distribution of rainfall and the annual temperature range at 100-year timeslices across the record. This model is corroborated by intervals of monthly-resolved laser ablation trace element data, providing glimpses into past Iberian seasonality shifts. The most salient features of the modelled results include extremely dry Younger Dryas winters (~12.9e11.6 ka BP) and several intervals during the midHolocene with almost no summer rainfall (e.g., at 4.2 and 9.0 ka BP). By 1.6 ka BP, a near-modern rainfall seasonality was established. According to the modelling results, seasonal rainfall and temperature distribution variability can account for 95% of the record. The model presented here provides a new tool for extracting critical missing seasonality information from stalagmite d18O records. Intervals where the model does not converge may represent transient climate anomalies with unusual origins that warrant further investigation. © 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
The Iberian Peninsula is located in a climatologically important location at the western periphery of Europe. The region's climate is sensitive to the phase of the North Atlantic Oscillation (NAO) (Rodrigo et al., 2000;Trigo et al., 2004) as well as oceanic heat advection to the North Atlantic due to Atlantic Meridional Overturning Circulation (AMOC) variability (Mary et al., 2017). The Iberian Peninsula is affected by both Atlantic and Mediterranean influences (Gimeno et al., 2010), and therefore the entire peninsula does not experience identical climate shifts at the same time; instead climate shifts are spatiotemporally variable. However, the proximity of northern Iberia to the climatologically-important oceanic subpolar and subtropical gyres (that transfer heat and salt toward the Nordic seas) (Morley et al., 2011;P erez-Brunius et al., 2004) and the Azores High (that partially controls the position and strength of westerly winds and the North Atlantic Oscillation) (Baker et al., 2015;Olsen et al., 2012;Trouet et al., 2009;Walczak et al., 2015) means that climate records from the region are particularly sensitive to major modes of atmospheric and oceanic circulation that also affect the rest of Europe and the entire Northern Hemisphere. The coastal region north of the Cantabrian Mountains exhibits particularly robust and stationary relationships between the NAO state and winter rainfall amount, at least over the most recent 100-year period for which climate reanalysis data are available (Comas-Bru and McDermott, 2014). AMOC variability, linked via teleconnections to the Intertropical Convergence Zone (ITCZ), also affects the region's climate (e.g., Pohlmann et al., 2006;Souza and Cavalcanti, 2009). North Iberia is therefore an ideal location to study both temporal shifts in AMOC as well as NAOdriven changes in temperature and precipitation.
However, this same sensitivity to several key European climate forcings also complicates the interpretation of north Iberian climate proxy records. This issue is compounded by the fact that the controls on the climate signal within any one proxy are temporally variable. A number of high-quality climate records exist that use terrestrial archives distributed across northern Iberia, and these have helped shed light on this issue. In particular, the region's extensive karst has permitted the development of a series of excellent stalagmite-based proxy records of temperature (Martin-Chivelet et al., 2011), rainfall (Moreno et al., 2017;Railsback et al., 2011;Smith et al., 2016), and the d 18 O composition of North Atlantic surface water (Domínguez-Villar et al., 2009). Recently, a number of publications have highlighted that the complex and sometimes contradictory pattern of climate proxy results observed could stem from shifting seasonality in rainfall (Morellon et al., 2009;Moreno et al., 2010Moreno et al., , 2017Walczak et al., 2015), potentially due to meridional displacement of the Azores High which broadly tracked insolation and mean North Atlantic climate state. Specifically, research has highlighted the possibility that climate change on the Iberian Peninsula is the result of time-transgressive seasonality changes across the region (Moreno et al., 2017). Different cave sites and stalagmite samples are variably sensitive to local environmental factors, and consequently the absence of a single continuous record from a single stalagmite covering the entire Holocene and into the last glacial complicates efforts to isolate the seasonality signal inherent within north Iberian climate.
Here we present new decadal-scale oxygen isotope data for a Useries dated stalagmite from La Garma Cave in north Iberia covering the entire Holocene. These new data, combined with previously published data from the same stalagmite over the interval from 10.5 to 13.5 ka BP, provide a rare continuous record of Iberian climate since the Bølling-Allerød interstadial from a single stalagmite sample. Consequently, the new isotope dataset provides the opportunity for replication, continuation, and aggregation of previously published shorter but otherwise high-quality Younger Dryas and Holocene speleothem records from northern Spain (Domínguez-Villar et al., 2008Dominguez-Villar et al., 2017;Martin-Chivelet et al., 2011;Moreno et al., 2010Moreno et al., , 2017Railsback et al., 2011;Rossi et al., 2018;Stoll et al., 2009). The oxygen isotope results are clearly not explicable by appealing to a single control (i.e., exclusively temperature, rainfall amount, etc.). We therefore explore the potential role of seasonality shifts in rainfall on the new d 18 O record by developing a new modelling technique that uses the oxygen isotope data to estimate rainfall and temperature seasonality shifts across the Holocene and into the late Pleistocene. This model is corroborated by intervals of laser ablation trace element data obtained at a monthly resolution, which provide direct glimpses into past seasonality shifts. The new data and model presented here are discussed in both site-specific and regional contexts, and, in combination with other published records from the Atlantic margin of Europe, offer new insights into Iberian climate evolution from the late Pleistocene. Specifically, modelled seasonal temperature and rainfall shifts provide a framework for the interpretation of other proxy records.

Site description
La Garma Cave (43 25 0 N, 3 40 0 W) is developed within lower Cretaceous limestones on several levels within a 187 m high hill situated 11 km ESE of Santander, northern Spain (Fig. 1a). The cave's entrance is located~5 km inland from the Bay of Biscay at 85 meter above sea level (m.a.s.l.). This cave site is well studied, and a detailed site description is available in previous publications (Arias, 2009;Arias et al., 2001Arias et al., , 2011Arias and Ontañ on, 2012;Baldini et al., 2015). The current vegetation above the cave consists of dense C 3 vegetation, including kermes oak, hazel, bay, and eucalyptus (Rudzka-Phillips et al., 2013). Soil thickness varies, but is typically about 1 m. Meteorological data for the Santander Airport Global Network of Isotopes in Precipitation (GNIP) station (43 29 0 28 N; 3 48 W; 52 m.a.s.l.), located 13 km west of La Garma Cave, indicate a mean annual air temperature of 14.78 C and a mean total annual precipitation of 870 mm from 2000 to 2010, and the K€ oppeneGeiger climate classification is 'warm temperate' (Comas-Bru and McDermott, 2014). The Santander Airport GNIP station also provides over a decade of monthly precipitation d 18 O (d 18 O p ) data. An analysis of longer rainfall records (from 1912 C.E.) from several local meteorological stations (Santander/Parayas, Villaverde de Pontones, and Santander Ojaiz) suggests that the GNIP data were measured during a slightly drier than average interval of the last 100 years, with mean rainfall values from these stations suggesting a long-term (1912e2011 AD) mean annual rainfall of~1228 mm compared with a mean annual rainfall of 1109 mm during the same interval as the GNIP data (2000e2010 AD). Here we use the GNIP data because it is coupled with rainfall  (Baker et al., 2019), which for our site is predominantly a consequence of temperature and moisture source.

Sample GAR-01 description and preparation
Stalagmite GAR-01 consists entirely of coarsely crystalline calcite and was collected in two parts on separate occasions in 2004 from the cave's Lower Gallery. The Holocene portion of GAR-01 (GAR-01A) was found lying on the cave floor adjacent to the actively growing portion of the same stalagmite (GAR-01B), which was collected. Sectioning revealed that recent calcite was deposited  (Mary et al., 2017); Shackleton Site (Ausin et al., 2019) non-conformably on top of older growth within GAR-01B (Fig. 2). Subsequent U-Th dating revealed that GAR-01A was broken from its growth position during the Middle Ages, and represents the middle interval of the adjacent, actively growing stalagmite GAR-01B that was still in situ. Archaeological analysis and 18 AMS radiocarbon dates of charcoal (n ¼ 11) and human bone (n ¼ 5) samples collected from the Lower Gallery of La Garma Cave suggest humans visited the location at least twice within the period 688e754 C.E. (Ontañ on et al., 2018), consistent with stalagmite U-Th ages for the break. The cave contains numerous stalagmites that were broken during the Middle Ages for unknown reasons, many of which have actively re-grown since the breakage, resulting in post-Middle Ages calcite deposition atop the broken bases. The physically separated but intrinsically linked samples GAR-01 A and B are collectively referred to as stalagmite GAR-01 in the text, unless specifically referring to either of the two physically discrete samples.
Stalagmite GAR-01 grew continuously from~14.0 ka to the date of collection in 2004 and represents 800 mm of total growth with a 57 mm yr À1 mean extension rate. Both GAR-01A and B were sectioned, polished, cleaned, and sampled using a conventional dental drill at a resolution of~38 years per sample. The results of high-resolution (sub-decadal) laser ablation analysis for C and O isotopes over the Younger Dryas interval was reported previously .

GAR-01 stable isotope analysis
Conventionally drilled powder samples (drilled at a 2.2 mm spatial resolution) were analysed using a GV Instruments Multiflow-Isoprime systems at Royal Holloway University, London (RHUL). Herein, the GAR-01 d 18 O record will be discussed in detail. Previous research has shown stalagmite d 13 C to be sensitive to  (in ka BP 1950). U-Series powders were drilled at regular intervals along the length of GAR-01. GAR-01 A was the first section of stalagmite GAR-01 to be collected and dated. GAR-01 A was broken off by Middle Ages visitors to the cave and discovered on the ground adjacent to the in situ portion of GAR-01 (GAR-01 B) in 2003. Thus stalagmite GAR-01 A constitutes the Holocene and GAR-01 B constitutes a record of the pre-Holocene with Middle Ages break and post-Middle Ages growth to present-day. The lower portion of GAR-01 B was the subject of a previous study that yielded a high resolution reconstruction of the Younger Dryas Event . GAR-01 (GAR-01 A and GAR-01 B) grew continuously from 13.864 ka BP to the date of collection in 2004. Some brief breaks in the record exist where a rock saw was used to slab the stalagmite for LA-ICPMS analysis. land-use changes at the surface (e.g., Baldini et al., 2005). Archaeological investigations at La Garma Cave over several decades by P. Arias and colleagues suggests that human activities (e.g., deforestation, fortifications, etc.) have intermittently modified the surface environment through the Holocene. Due to the potential for these activities to overprint the d 13 C climate signal, an in-depth discussion of the GAR-01 d 13 C record is considered outside the scope of the current study but will form the basis of a future paper on anthropogenic activity at the site.

LA-ICPMS analysis
A 193-nm wavelength excimer LA-ICPMS system at RHUL (RESOlution M-50 prototype coupled to an Agilent 7500ce quadrupole ICP-MS) featuring a two-volume Laurin LA cell was used to produce a high-resolution Mg dataset across select intervals of the stalagmite. Trace element data were obtained as continuous profiles using a rectangular spot (285 Â 12 mm) at a 10 mm s À1 sample movement and a 15 Hz laser repetition rate, producing an effective spatial resolution of~15 mm, which is approximately monthly-scale for fast growing intervals of GAR-01 during the mid-Holocene. 43 Ca concentrations were used as an internal standard and NIST612 as external standard (Müller et al., 2009).

U-Th dating and age model development
Twenty-four powder samples were drilled from within distinct growth layers along the central growth axis of GAR-01 using a handheld drill and a tungsten carbide drill bit for uranium series dating. Chemical separation and purification of uranium and thorium and their subsequent isotopic analysis using a Thermo-Finnigan Neptune multicollector inductively coupled mass spectrometer (MC-ICP-MS) at the University of Bristol followed the techniques and procedures outlined by Hoffmann et al. (2007). Measured 238 U concentrations in stalagmite GAR-01 ranged between 71.2 and 118.7 ng g À1 , and the measured 230 Th/ 232 Th activity ratio varied between 22.1 and 2089.9. All ages were calculated using the decay constants reported in Cheng et al. (2000) and corrected for detrital contamination assuming a bulk earth ( 238 U/ 232 Th) activity ratio of 0.746 (Hans Wedepohl, 1995) and ( 230 Th/ 238 Th) ¼ ( 234 U/ 238 U) (Table S1).
All 24 U-Th dates from GAR-01 are in stratigraphic order (Table S1) and a distance-age model was generated using the COPRA (Constructing Proxy Record from Age models) algorithm . For the previously published Younger Dryas record, we used the StalAge algorithm to calculate the age model , but here considering the proxy record with translated age model uncertainties facilitates comparisons to the large number of shorter existing records, and COPRA has this capability. Although both the COPRA and StalAge algorithms construct age models using Monte Carlo simulations to interpolate between U-Th ages and account for potential outliers, age inversions, and potential hiatuses, StalAge does not yet translate the age model uncertainties to the proxy record. Additionally, the StalAge mean age model for the full GAR-01 record contained short intervals where the StalAge model produced unrealistically high growth rates ( Supplementary Fig. S2). For these reasons we elected to use the COPRA algorithm in the present study. However, the two age models agree very well overall, particularly over the Younger Dryas interval (Supplementary Fig. S2). Based on 24 230 Th dates, the 2004 C.E. date of collection, and the COPRA output, stalagmite GAR-01 grew continuously from 13.861 ka until it was collected in 2004 C.E. (Fig. 3 and Supplementary Table S1). Based on the age model, the conventional drill d 18 O dataset discussed herein has a mean temporal resolution of 37.9 years.

Cave monitoring
Air temperatures were monitored at four different locations within La Garma Cave from July 2006 to July 2007 using TinyTag temperature loggers (Jackson, 2009), including at the location from which stalagmite GAR-01 was collected. Monitored cave air temperatures at all of the La Garma Cave monitoring sites are 1.38e2.70 C lower than those of the surface meteorological measurements. The elevation of the monitored cave passage (59 m.a.s.l.) is close to that of the meteorological station (52 m.a.s.l.) and therefore altitude cannot account for this difference. Three of the monitored sites in the cave recorded air temperatures with extremely small variations through the year (annual range < 0.22 C). However, air temperature at the GAR-01 site has a higher annual range of 0.69 C (from a maximum monthly mean value of 12.37 C in November and a minimum of 11.68 C in February), and a long-term mean temperature value of 12.08 C. This is 2.70 C colder than the mean outside temperature, and reflects cave ventilation systematics and a winter cold air trap, common to temperate cave sites (Bourges et al., 2001;James et al., 2015;Sherwin and Baldini, 2011), that disproportionately affects the GAR-01 site.
As previously reported , drip samples were obtained during monthly site visits to La Garma cave between September 2004 and October 2005. Monthly-integrated dripwater samples were collected from: i) the feeder drip site to stalagmite GAR-01 (GDW-1) and ii) a drip site (GDW-2) located 3 m from drip site GDW-1. Dripwater samples were collected using 1 L graduated collection vessels to allow the mean drip rate to be calculated for each period of collection. More recently (June 2012 to November 2013), and to augment the monthly-scale drip rate data, a Sta-lagmate® drip logger was deployed at the GAR-01 drip site to monitor drip rate at 30-min intervals over 18 months. All dripwater samples were analysed for oxygen and hydrogen isotope compositions. All drip water isotope analyses were performed at the Nevada Stable Isotope Laboratory, University of Nevada-Reno (UNR) using a Micromass Aquaprep device interfaced to a Micromass dual inlet IsoPrime stable isotope ratio mass spectrometer and the Epstein and Mayeda (1953)  Stalagmite GAR-01 was fed by a seasonal drip on the basis of a maximum discharge rate of 4.09 Â 10 À7 L s À1 and a coefficient of variation of 55.5 (i.e., calculated as the absolute value of the standard deviation divided by the mean drip rate multiplied by 100) following the Smart and Friederich (1987) classification scheme. Stalagmate® drip logger data (collected between June 2012 and November 2013) revealed that the GAR-01 drip site (GDW-1) drips continuously (i.e., never drying completely) and exhibits long-term seasonal trends with daily drip rates of 131 drips per day in June 2012 decreasing to a minimum of 75 drips per day in winter 2012.
The mean GDW-1 dripwater d 18 O value (À5.61‰, V-SMOW) combined with the d 18 O value of the calcite that forms the most recent deposition on stalagmite GAR-01 (À3.93‰, V-PDB; the last drill analysis point, centred around 1976 C.E.) were used to assess the extent to which GAR-01 calcite was deposited in oxygen isotopic equilibrium with its dripwater. Under perfect equilibrium conditions, the temperature calculated on the basis of equilibrium water-calcite fractionation factors should approximate the observed modern cave air temperature. However, considerable debate exists in the literature regarding the most appropriate water-calcite oxygen isotope fractionation factor for speleothem calcites (e.g., Demeny et al., 2010;Fairchild and Baker, 2012;Johnston et al., 2013;McDermott, 2004;McDermott et al., 2011;McDermott et al., 2006;McDermott et al., 2005;Tremaine et al., 2011), and even regarding the extent to which the equilibrium concept is applicable to stalagmites (Da€ eron et al., 2019). Use of the Kim and O'Neill (Kim and O'Neill, 1997;Kim et al., 2007) equation yields a temperature of only 6.3 C, much lower than the measured mean cave air temperature at the GAR-01 site (c. 12 C), reflecting the well-known tendency for this equation to yield apparent cave temperatures that are much too low (e.g., McDermott et al., 2005). By contrast, the empirical equation of Tremaine et al. (2011), based on 'farmed' speleothems at multiple locations within Hollow Ridge Cave (Florida), indicates a value of 10.6 C, a little lower than the measured cave air temperature. The equations of Coplen (2007) and Demeny et al. (2010) yield values that are higher (13.29 C)  . The mean growth rate over the Holocene was 59 mm yr À1 . Intervals with particularly rapid growth rates are observed during the late Northgrippian stage (between 4.7 and 4.2 ka BP: 156 mm yr À1 ) and during the Greenlandian stage (between 10.4 and 8.7 ka BP: 149 mm yr À1 ). The intervals of slowest growth occur during the Younger Dryas (discussed in detail in Baldini et al. (2015)) and during the mid-Meghalayan stage (between 3.2 and 1.8 ka BP: 14.6 mm yr À1 ). and lower (9.0 C) than the observed cave air temperature, respectively. Clearly it is not possible to be definitive about the extent to which cave carbonates are precipitated in oxygen isotopic equilibrium and our perception of equilibrium depends on the somewhat arbitrary choices of fractionation factors from the literature. Regardless, the fractionation factor based on 'cave-farmed' speleothems (Tremaine et al., 2011) is found to most closely approximate the observed cave air temperature at La Garma, suggesting that this is most appropriate for GAR-01 calcite deposition. temperature (both regional external air temperature and in-cave air temperature), rainfall amount (the 'amount effect'), moisture source isotopic composition, rainfall seasonality, and moisture mass trajectory, with temperature and trajectory dominating the signal in northern Iberia. Correspondingly, variability across the and November 2005 (blue circles) and La Garma cave drip water from drip sites, GDW-1 (feeder drip to stalagmite GAR-01; light green circles) and GDW-2 (a drip site 3 m away from GDW-1; green crosses) sampled during the study period. B, C) Dripwater isotope data for GDW-1 (light green circles) and GDW-2 (green crosses) compared to contemporaneous estimates of drip rate (based on dividing the volume of drip water collected by the period of collection) through time. D) Variability in effective precipitation (an estimate of the net hydrological balance) and Santander rainwater d 18 O during the study period. E) Variability in monthly actual evapotranspiration and Santander precipitation amount during the study period. Actual evapotranspiration and effective precipitation were calculated using the modified Thornthwaite (Thornthwaite, 1955) and Grindley (Grindley, 1969) methods. Horizontal dashed grey lines indicate the amount-weighted mean (WM) of drip water d 18 O for GDW-1 (À5.73‰) and GDW-2 (À5. Younger Dryas (YD) interval in GAR-01 was previously interpreted as primarily reflecting external temperature (i.e., low regional temperatures driving meteoric precipitation d 18 O and consequently calcite d 18 O to lower values) with moisture source region and trajectory shifts playing a secondary role . According to Baldini et al. (2015), a temperature decrease of 6e9 C during the YDE compared to Bølling-Allerød (B-A) temperatures could explain the 3.1‰ lowering of d 18 O during the event. Elevated ocean water d 18 O and lower in-cave temperatures characteristic of that time interval would have forced the GAR-01 d 18 O record higher, so cannot explain the observed decrease, although both may have offset the observed decrease somewhat. Similarly, sea level rose~40 m across the YD from 14 to 11 ka BP (Grant et al., 2014), and a more distal shoreline similarly cannot explain the lower values. However, the new 14 ka-long GAR-01 record presented here reveals d 18 O values during several intervals within the Holocene (an epoch characterised by warm conditions in northern Iberia) that are similar to those observed during the YD (Mary et al., 2017). For example, in the decadally-resolved GAR-01 d 18 O dataset, the d 18 O at 8.974 ka BP is À5.08‰, compared to a value of À5.01‰ at 12.29 ka BP at the height of the YD; similarly the d 18 O at 4.709 ka BP is À5.00‰ (Fig. 5). The Iberian Margin sea surface temperature (SST) record from off the SW coast of Iberia (39 11 0 N, 10 0 0 W) suggests a mean early Holocene SST of~19 C (Bard, , 2003Pailler and Bard, 2002a, b) compared with a YD low of~13 C. A SST record from the Bay of Biscay (Mary et al., 2017) suggests that north Iberian temperatures at~4.7 ka BP were~14.9 C, warmer than the modern values from the same core of 14.7 C, and almost certainly warmer than YD values (the Bay of Biscay SST record does not extend to the YD). This is also consistent with other regional records (such as Greenland ice core temperature proxy records) has a columnar and crystalline structure with an absence of visible laminae (Fig. 2), and iv) the~6e7.5 month residence time (Fig. 4) combined with a lack of dissolution cups in the stalagmite morphology suggests that the drip at no time becomes undersaturated with respect to calcite.

GAR-01 stable isotope data across the Holocene
The GAR-01 conventionally drilled d 18 O data between 13.345 ka BP and 2004 C.E. ranged from À5.08 to À3.34‰ (mean ¼ À4.37‰, 5b). For clarity, the Holocene GAR-01 record is subdivided into the following three Holocene stages as recently defined by the International Commission on Stratigraphy (ICS) (Cohen et al., 2013;Walker et al., 2012): i) Greenlandian (11.7 to 8.2 ka BP), ii) Northgrippian (8.2 to 4.2 ka BP), and iii) Meghalayan (4.2 ka BP to Present). During the Greenlandian stage (11.7 to 8.2 ka BP), the most conspicuous feature is a decreasing trend in the GAR-01 d 18 O data from À3.84‰ at 11.7 ka BP to À4.99‰ at 8.2 ka BP ( Fig. 5). The most pronounced negative excursion during this interval occurs at~9 ka BP and may reflect the north Iberian climatic response to the previously documented '9.2 ka event'. Detected in numerous climate proxy records across the Northern Hemisphere (Fleitmann et al., 2008;Genty et al., 2006;Masson-Delmotte et al., 2005;Rasmussen et al., 2007), the '9.2 ka event' may reflect the effects of a meltwater pulse (MWP) (Fleitmann et al., 2008). The event is also observed in a stalagmite record from Dongge Cave in China (Dykoski et al., 2005), suggesting that it also influenced the East Asian Summer Monsoon. Interestingly, the expression of the 9.2 ka event in GAR-01 is similar in magnitude to the GAR-01 YD d 18 O anomaly  with d 18 O values of À5.08 and À5.00‰, respectively (Fig. 5). Previous research suggests that conditions in northern Spain were likely considerably warmer and wetter at~9 ka BP (Morell on et al., 2018) relative to the YD and this is consistent with relatively rapid GAR-01 growth rates observed during the early Holocene (Fig. 3) (Bond et al., 1997;Bond et al., 1999); however, the isotope response within GAR-01 is relatively subdued. Our data suggests that the shift to drier conditions at 4.2 ka BP at this North Atlantic maritime site was muted compared to the pronounced aridification that occurred across other drier areas (Carolin et al., 2019;Isola et al., 2019;Staubwasser et al., 2003). This interpretation is further supported by the relatively rapid GAR-01 growth rate across the latter 500 years of the Northgrippian stage (Fig. 3). The Meghalayan stage (4.2 ka BP to 2004 C.E.) in stalagmite GAR-01 is characterised by a 0.45‰ negative shift between 4.0 ka BP and 3.7 ka BP followed by a gradual increasing trend to À3.96‰ (at 1.8 ka BP  (Fig. 5) exhibit a pronounced negative anomaly (À0.67 and À1.13‰, respectively) that is not associated with any known climatic event, nor is it replicated in regional records (Fig. 6). Thus, we tentatively ascribe the anomaly to site-specific factors (e.g., a localised climate shift or hydrological event) rather than a regional scale climate event.  Martin-Chivelet et al., 2011;Moreno et al., 2010Moreno et al., , 2017Railsback et al., 2011;Rossi et al., 2018;Smith et al., 2016;Stoll et al., 2009) provide valuable overlapping data coverage and a degree of replication for sections of the GAR-01 Holocene time-series (Fig. 6). In Fig. 6 between GAR-01 and regional records (d 18 O STAL ) largely reflects differences in cave altitude and distance from the coast (Fig. 1b). Stalagmites from all sites plot along statistically significant negative linear regression lines when the mean d 18 O offset is compared to cave distance from the coast (y ¼ À29.75x e 0.60; r 2 ¼ 0.96, p < 0.001) and altitude (y ¼ À359.80 x þ 92.70; r 2 ¼ 0.91, p < 0.001) except the ESP03 record from Cova da Arcoia that has mean d 18 O offset (À0.45‰) relative to GAR-01 that is higher than expected given the cave's high altitude and distance from the coast (Fig. 1b). The slope of the isotope ratio versus altitude line is similar to that calculated by de Oliveira and da Silva Lima (2010)   and distance from the coast of Molinos and Ejulve caves relative to La Garma Cave (Fig. 1b).
To determine the relationship between regional and GAR-01 d 18 O records, all datasets were linearly interpolated using the Regular Interpolation (linear) function within the statistics software PAST v. 3.x (Hammer et al., 2001) and temporally overlapping intervals were compared using Spearman's rank correlation analysis.
On decadal timescales regional d 18 O datasets are significantly but weakly correlated with GAR-01 d 18 O, despite the strong visual match (Fig. 7). The lack of stronger correlations is probably due to chronological uncertainties and a generally low signal-to-noise ratio characteristic of Holocene climate (Figs. 6 and 7). For instance, GAR-01 is negatively correlated with Asiul stalagmite ASM (Spearman's rank correlation coefficient (r s ) ¼ -0.20, p < 0.001) whilst ASR from the same cave exhibits a significant positive correlation (r s ¼ 0.24, p < 0.001) with GAR-01 (Fig. 7). Additionally, stalagmites ESP03 and GAR-01 are significantly positively correlated (r s ¼ 0.21, p < 0.001) over the pre-4 ka BP portions of these records but the large dating uncertainties of the ESP03 record over the Meghalayan (i.e., the Late Holocene) prevent statistical comparison over that interval (Fig. 7). Besides dating uncertainties, other factors potentially responsible for the weak correlations observed between GAR-01 d 18 O and regional records at decadal timescales include site-specific differences (e.g., varying hydrologies lending certain samples more or less sensitive to transient climate events ) and possible aliasing effects (e.g., differences in the temporal resolution of sampling; Baldini et al., 2008). Although GAR-01 and Katie Cave's LV5 record exhibit similar variability (r s ¼ 0.20, p < 0.01) between 8.545 and 6.905 ka BP (with both detecting a two-pronged 8.2 ka BP event), the overall correlation of Kaite cave stalagmites (LV3 e 6) with GAR-01 is quite low (r s ¼ 0.09, p < 0.01). This is potentially due to chronological uncertainties, or alternatively to Kaite Cave's location on the southern (leeward) slopes of the Cantabrian Mountains that form a hydrological/meteorological divide between northern and central Spain (Fig. 1a).

Model set-up
As discussed previously, the variability within the GAR-01 d 18 O record is not reasonably attributable exclusively to mean annual temperature shifts. In particular, low values of~-5‰ found at 12.29 ka BP (the height of the cold YD), at 8.974 ka BP (regional warmth), and at 4.709 ka BP (regional warmth) cannot have the same origin. However, it is conceivable that seasonality changes in both rainfall and temperature may have contributed to these common low values. To investigate this further, we have developed a simple model that uses a Monte Carlo approach to estimate the seasonal distribution of rainfall as well as the annual temperature range. The model assumes that annual mean regional outside temperature is represented by the Iberian Margin temperature record (Bard, , 2003) adjusted for regional temperature differences (Table S2). The correlation between modern surface air temperatures between the nearest meteorological stations to both sites (La Garma Cave and the Iberian Margin core site) is very high (Fig. S3). The Mary et al. (2017) Bay of Biscay SST record is more proximal and higher resolution than the Iberian Margin record (Fig. 1a), but does not extend into the YD, and is therefore compared with the output of the model rather than as an input. Modern temperatures at the La Garma Cave site are~3 C cooler than at the Iberian Margin core location, consequently modelled MAT outside the cave is corrected by 3 C compared with the Iberian Margin SST record. The annual in-cave temperature is set as 2.7 C colder than the regional mean annual outside temperature. Whereas at many sites in-cave temperature reflects the outside MAT, at La Garma Cave monitoring reveals the presence of a 'cold trap' that exaggerates the importance of winter outside air temperatures (see Section 3.5). Monthly modelled temperatures are linked to the outside mean annual temperature via modern empirical relationships between annual mean, maximum, and minimum temperature values observed at the Santander GNIP station. The modern annual temperature range (e.g., seasonality) increases only slightly with decreasing temperature, and is controlled by the positive relationship between the mean annual air temperature and the minimum and maximum temperature, respectively, in which the minimum temperature is more sensitive to changes in the mean annual temperature (Table S2). The model permits a randomly generated maximum increase or decrease in annual temperature range of 1 C; in other words, minimum and maximum temperatures derived from the Iberian Margin SST record are allowed to randomly vary by between À0.5 and 0.5 C each, and consequently the maximum range introduced into the modelled temperature output range above and beyond that observed in the Iberian Margin SST record is 1 C. Consequently, the modelled MAT also has a maximum randomness of 1 C. Because the maximum randomness in temperature range introduced by the model is limited to 1 C, the model returns whether the seasonal temperature range increases or decreases, but may underestimate the actual increase in seasonality. The model assumes that the same relationship between meteoric precipitation d 18 O and temperature observed at the nearest GNIP station (Santander, 13 km to the west of the site) holds true through the entire GAR-01 record (see Table S2 for a detailed list of parameters and steps used in the model). Next, the model creates a seasonal rainfall distribution of random polarity and amplitude. In other words, the model produces a randomly generated seasonal cycle in rainfall that could range from an extremely wet winter but arid summer, to an extremely arid winter but wet summer, to extremely muted seasonality (i.e., an even distribution of rainfall through the year). The model then produces an amount weighted annual mean meteoritic

Assumptions within the model
The model contains several assumptions. For example, the total annual rainfall amount is held constant, whereas rainfall amount undoubtedly varied interannually. This assumption is inconsequential however because rainfall seasonality drives dripwater d 18 O for well-buffered sites rather than the total annual rainfall; regardless, the model cannot return annual rainfall totals. At this site, the predominant control on rainfall d 18 O is temperature (i.e., the season in which the rain occurs) rather than the rainfall amount. This is supported by a strong relationship between monthly rainfall d 18 O and temperature at the Santander GNIP site (r s ¼ 0.56; p < 0.001; n ¼ 130), however, varying moisture source region, trajectory, and moisture source d 18 O undoubtedly affect rainwater d 18 O as well. The model assumes that rainwater d 18 O is driven entirely by local temperature variability. Vegetation above La Garma Cave is also held constant throughout the model, whereas The GAR-01 record is continuous, but only data overlapping with the other records is shown for clarity, resulting in apparent gaps in the GAR-01 record. Note the ASM y-axis is inverted because the r value is negative (potentially due to chronological uncertainty). In panel c, the GAR-01 record is normalised to zero to facilitate comparison with the normalised Asiul Cave record. in reality climatologically-and anthropogenically-driven vegetation shifts above the cave probably occurred and affected seasonal infiltration amounts (Baldini et al., 2005). Although local pollen records provide information regarding regional vegetation type through the Holocene, a localised reconstruction of vegetation density immediately above the cave has yet to be attempted, thus, this parameter is currently unquantifiable. The model includes options to parameterise soil evapotranspiration using either the Thornthwaite equation or Hamon's equation. However, the model results are optimised only when evapotranspiration is set to zero importance. 'Switching off' evapotranspiration in the model is justified based on recent empirical soil-water d 18 O results from La Garma Cave where evaporation was found to be negligible (Comas-Bru and McDermott, 2015). However, transpiration may play a role insofar as it can cycle moisture from the soil directly back to the atmosphere and if this is seasonal (likely for deciduous vegetation) it could impact the seasonal moisture available for infiltration to the cave. Additionally, over relatively long periods of the Holocene, vegetation changes (e.g. density, moisture use efficiency, etc.) may also have an effect. Future versions of the model may include the capability of incorporating a vegetation model pending the acquisition of new vegetation density data for the site (e.g., from local pollen reconstructions or organic markers in speleothems).
The model also assumes that the kinetic effects that are inherent to the Tremaine water-calcite fractionation factor are constant yearround; because no research has yet quantified this variability at this particular site, we are not able to incorporate this into the model. However, future versions could incorporate assumptions regarding the seasonality of disequilibrium effects (Deininger et al., 2012;Deininger and Scholz, 2019;Mühlinghaus et al., 2007) on predicted responses of the cave environment to the external climate signal, ideally by incorporating new models such as IsoCave (Guo and Zhou, 2019) or ISOLUTION (Deininger and Scholz, 2019). Research suggests that growth rate affects oxygen isotope ratios in calcite (Hansen et al., 2019;Stoll et al., 2015), although we note that these interpretations have been challenged (Dreybrodt, 2016;Dreybrodt, 2019). This first version of the model does not calculate seasonal growth rates, so it is not possible to correct for possible growth rate effects yet, although this is planned for a future version of the model. However, the observed relationship between growth rate and d 18 O may reflect a common control on both rather than a cause and effect relationship between growth rate and d 18 O (Dietzel et al., 2009;Fohlmeister et al., 2018;Gabitov et al., 2012) However, equally informative are the intervals when the model cannot successfully converge on a value, indicating intervals when seasonality shifts are insufficient on their own to produce the stalagmite d 18 O values. In these few instances, another factor is necessarily implicated, highlighting intervals with considerably different boundary conditions that require a different interpretation. In theory the model could produce results on any timescale longer than annual. Here, we use 100-year long timeslices obtained by interpolating the decadally-resolved d 18 O data at regular 100year intervals; the process repeated for each timeslice.

The temperature model
The modelled monthly temperature output for the most recent timeslice matches observed monthly distribution very well, suggesting that the model is robust. The YD interval of the GAR-01 record is discussed qualitatively in Baldini et al. (2015), but here we provide quantitative estimates of seasonal temperature and rainfall (section 5.3.2. below). The YD interval within the modelled temperature data is characterised by winter temperatures that arẽ 7 C lower than modern winter temperatures, and summer temperatures that are~4 C lower than modern summer temperatures (Fig. 8). However, our approach is limited to a maximum of 1 C of seasonal temperature change above and beyond that suggested by the Iberian Margin SST record (Bard, , 2003, so the actual temperature range may exceed this estimate. Here we do not attempt to reconstruct absolute temperatures, but instead highlight the seasonal differences in temperature and rainfall that could account for the observed d 18 O signal. Despite uncertainty regarding the absolute temperature range, this is consistent with a number of YD reconstructions suggesting that the event was characterised by considerably lower wintertime temperatures (Denton et al., 2005;Lie and Paasche, 2006), possibly induced by an AMOC slowdown (Lynch-Stieglitz, 2017). The recovery out of the event was rapid, and according to the GAR-01 modelling, by 9.4 ka BP monthly temperatures were slightly warmer than present day (by~0.5 C) (Fig. 8). The rapid amelioration in temperature at the site likely reflects the well-documented rapid deglaciation following the YD and the maximum in 60 N summer insolation centred on~10.5 ka BP (Fig. 5) (Berger and Loutre, 1991). From the early Holocene to the date of collection, modelled temperatures gradually cool with modern winter and summer temperatures~0.5 C and~0.3 C cooler than early Holocene values. The gradual Holocene cooling, however, is punctuated by a series of warming and cooling events. Most notably, the 8.2 ka event is clearly expressed in the GAR-01 d 18 O record, and modelled temperatures reach a local minimum during the 8.2 ka timeslice, with winter temperatures of 10.0 C (compared with 11.7 C during the early Holocene). This suggests that the low d 18 O values observed in the GAR-01 d 18 O dataset resulted at least partially from cooler temperatures and increased seasonality (see section 5.3.2. below), consistent with previous results from other Atlantic margin sites Daley et al., 2016;Domínguez-Villar et al., 2009). This is the lowest winter temperature modelled for the early Holocene, and in fact only at 4.3 ka BP do winter temperatures drop below the 8.2 ka values, reaching as low as 9.5 C, potentially reflecting the welldocumented '4.2 ka event' (Fig. 8). The expression of the '4.2 ka event' in the d 18 O data is rather muted, although it does appear in the Smith et al. (2016) record from nearby Asiul Cave. The event may have contributed to severe Middle East drought affecting civilisations at the time (Cullen et al., 2000;Hsiang et al., 2013). The expression of the event as one of the coolest timeslices of the Holocene, combined with its presence in other north Iberian records, suggests that the 4.2 ka event did affect the western margin of Europe as well as the Middle East. This is discussed in more depth below.
Modelled GAR-01 temperature variability is of a very similar amplitude as the Martin-Chivelet et al. (2011) temperature reconstruction over the last 4 ka (Fig. 8), with several key features duplicated. For example, both records suggest substantial centennial-scale temperature variability from 4 to 2.5 ka BP, and exhibit a pronounced warming event at~3 ka BP. A very pronounced warming event in the Mary et al. (2017) Bay of Biscay SST record at~2 ka BP is not reproduced either in the GAR-01 modelled temperature record or in the Martin-Chivelet et al. (2011) temperature reconstruction, suggesting that this temperature anomaly was largely restricted to the marine environment and only briefly affected terrestrial Iberian temperatures. Modelled GAR-01 temperatures show a maximum value at 1 ka BP, consistent with both the Martin-Chivelet et al. (2011) temperature reconstruction for northern Iberia (Fig. 8) and a comprehensive review of Iberian terrestrial climate proxy evidence (Moreno et al., 2012), and was probably linked to the Medieval Climate Anomaly (MCA). Over the most recent millennium, modelled GAR-01 temperatures suggest that minimum temperatures associated with the Little Ice Age (LIA) occurred from 1700 to the 1800 s C.E., rather than~300 years earlier as suggested by the Martin-Chivelet et al. (2011) temperature reconstruction. The modelled temperature values are however consistent with the most substantial advance in mountain glaciers across Iberia during the late LIA (Trueba et al., 2008), as well as extensive periglacial landforms dated to between 1700 and 1900 C.E. (Oliva et al., 2016(Oliva et al., , 2018. The modelled temperature output therefore appears to accurately reflect known temperature fluctuations during the last two millennia, capturing both LIA cooling and MCA warming. Overall, the excellent agreement between the modelled temperature output and existing temperature reconstructions over the last 4 ka BP strongly suggest that the modelling approach is reasonable, supporting interpretations over the older less-well constrained intervals of the record.

The rainfall model
The model run for the most recent timeslice (from 2000 to 1900 C.E., or À0.050 to 0.050 ka BP) converges on an annual distribution of rainfall which matches current observed seasonality (Figs. 9 and 10a), supporting the modelling approach. The model suggests that the seasonal rainfall pattern in northern Iberia had a modern polarity (i.e., rainier winters, drier summers) for 71% of the Holocene. Several intervals of reversed polarity (i.e., the opposite of present day) seasonal rainfall are also present (Fig. 9), and the model helps ground-truth interpretations based on the d 18 O datasets. One of the most notable intervals of reversed polarity seasonal rainfall occurs during the 12.80 ka BP timeslice, during the early YD (Fig. 10f). The modelled results suggest that peak rainfall at 12.80 ka BP occurred during the summer (~115 mm in August) whereas the driest month was November with~65 mm of rainfall (Fig. 10f). A very low GAR-01 d 18 O value exists at this date, which was interpreted by Baldini et al. (2015) as reflecting cold conditions consistent with an AMOC slowdown and increased sea ice across the North Atlantic. As discussed above (section 5.2.1.), the model supports this interpretation, but here we can add additional detail. If modern rainfall  (2002) (dashed white line). The offset between MAT at the La Garma site and Iberian Margin SST record site is corrected for using the most recent value in the Iberian Margin record, and this offset is applied throughout the rest of the study interval; this produces the slight apparent offset between the two records through the Holocene. The inset shows the modelled GAR-01 MAT and a temperature reconstruction for northern Iberia based on stalagmite d 13 C values (Martin-Chivelet et al., 2011), while the black crosses show a recent alkenone-derived SST record from offshore SW Iberia for comparison (SHAK06e5 K at 37 34N, 10 09W) (Ausin et al., 2019). Both the Martin-Chivelet and Ausin records are independent of the modelling, and have not been rescaled, so regional temperature differences remain. In contrast, by 9.4 ka BP temperatures were high (about 0.5 C higher than modern values), and the model suggests that summer rainfall was essentially non-existent ( Figs. 9 and 10d). An alternative is that summer rainfall did occur, but substantial amounts of summer evapotranspiration eliminated most summer recharge. However, this is not supported by the modelling results, where the incorporation of evapotranspiration reduces the ability of the model to converge onto a solution. This summer aridity suggests that a Mediterranean-like climate may have existed in northern Iberia at this time, consistent with a northward displaced Azores High that may have accompanied increasing insolation (Fig. 11). This is also consistent with the interpretations of Walczak et al. (2015) who argued that increased year-round moisture in southern Iberian was triggered by a northward displaced Azores High. The westerly storm track would have also shifted to the north, and indeed records suggest increased moisture delivery to Scandinavia (Bakke et al., 2009) and reduced moisture to southern France (Genty et al., 2006;Wirth et al., 2013) and northern Iberia (Aranbarri et al., 2014;Railsback et al., 2011) at this time. During the 8.2 ka event, winters appear wetter than at present day, but numerous other intervals also exist across the Holocene with similar amounts of winter rainfall. The model suggests that rainfall was strongly biased towards the summer over a 2 kyr period from 7 to~5 ka BP (Fig. 11). For example, at 5.9 ka BP, July rainfall was  and November (49.5 mm and 118.8 mm,respectively). In all cases, the total annual rainfall is held constant (consequently, the length of the bar for every timeslice is identical). The darker shading in both summer and winter indicates timeslices when rainfall totals are higher than current values.
153 mm compared to the modelled November total of just 31 mm (Fig. 10c). The end of this mid-Holocene interval of reversed polarity rainfall seasonality coincides with the establishment of Mediterranean conditions in southern Iberia at 5.3 ka BP, as well as the shift from a mostly positive NAO to a more negative NAO (Olsen et al., 2012). Both of these phenomena are likely due to southward displacement of the Azores High following the Holocene thermal maximum, consistent with elevated summer rainfall at 5.9 ka BP in northern Iberia slowly giving way to drier summers by~4.8 ka BP (Fig. 9). According to the rainfall model output (Figs. 9 and 10b  culmination of a summer drying trend starting at peak summer rainfall during the mid-Holocene; this biases the annual recharge towards low d 18 O winter rainfall. The well-documented '4.2 ka event' is linked to a reduction in Northern Hemisphere monsoonal system strength (Booth et al., 2005;Staubwasser et al., 2003) and is implicated in dramatic cultural change at several locations (deMenocal, 2001;Stanley et al., 2003), but particularly in the Middle East (Arz et al., 2017;Drysdale et al., 2006;Staubwasser et al., 2003). The event coincides with a peak in North Atlantic ice rafted debris possibly induced by solar variability (Bond et al., 2001). However, in our modelled rainfall seasonality dataset and in the d 18 O data, the 4.2 event summer aridity is the culmination of a longer summer drying trend starting after peak mid-Holocene summer rainfall at 5.9 ka BP (Fig. 9). This perspective is consistent with data from an Italian flowstone record also suggesting that the 4.2 ka event occurred near the end of a longer drying trend and that it may have a different origin than other Holocene drying events (Drysdale et al., 2006). Another interval of higher summer rainfall occurs from~2.8 to~1.8 ka BP, but rainfall occurs in the winter as well (although slightly reduced compared to modern values). From 1.6 ka BP to present day the rainfall distribution has a modern polarity, with the LIA characterised by enhanced seasonality with drier than modern summers. The rainfall and temperature distribution of the LIA is very similar to that observed in the model output during the 4.2 ka event, perhaps indicative of a similar origin.

LA-ICPMS analysis of the 4.2 ka event
Monthly-scale LA-ICPMS trace element data were obtained across the interval containing the 4.2 ka event, from 3.79 to 4.98 ka BP (Fig. 12). These data permit not only the reconstruction of climate with excellent detail, but also help corroborate the results of the seasonality model discussed above. Unfortunately determining the polarity of the seasonality is not possible, but the amplitude of any seasonality present is assessable.
Mg concentrations in GAR-01 were previously interpreted as reflecting seaspray contribution. Mg and Sr data over the 4.2 ka event in GAR-01 are consistent with this interpretation and an approximate marine aerosol contribution to drip water and, ultimately to GAR-01 calcite, of 2% (Fig. S4). Thus, aerosol-derived Mg concentrations decrease gradually from a peak at~5.0 ka BP until 4.2 ka BP, interrupted by several large Mg excursions (including at 4.9 ka BP) (Fig. 12a). This may reflect the decreasing influence of seaspray, following on from the interpretations of Baldini et al. (2015) for the YD interval within the same stalagmite. The mid-Holocene interval analysed containing the 4.2 ka event generally has considerably less Mg (545 ppm) than the mean value for the YD interval (843 ppm), potentially reflecting weaker winds and less seaspray affecting the coastal La Garma Cave site during the mid-Holocene than during the YD. Another possible factor is a shift in the direction of predominant winds due to a shifted Azores High, leading to a reduction of winds coming off the sea and consequently in seaspray. The model suggests pronounced summer aridity in northern Iberia between 4.2 and 4.5 ka BP, coincident with the lowest Mg values of the interval (<300 ppm). The mechanisms through which summer aridity might produce low Mg concentrations in GAR-01 include: i) reduced seaspray contributions due to a generally amenable climate and calm winds, even in winter, ii) very low Mg during winter due to high winter rainfall, or iii) a combination of i and ii.
During the YD interval, seasonality is generally not apparent in the Mg data (Fig. 12e), almost certainly reflecting low growth rates across this interval, where even the high-resolution LA-ICPMS data cannot tease out annual cycles. The growth rate during the YD interval is estimated at~25 mm year À1 and the LA-ICPMS data were obtained at a~15 mm effective resolution; consequently discerning an annual cycle is not possible. The mean Mg concentration values are considerably higher than the Holocene values, interpreted as reflecting more winter storms, stronger winds, and higher seaspray contributions to rainwater . A similar lack of seasonality is apparent at~4.8 ka BP (Fig. 12d), but the U-Th chronology suggests that the growth rate was sufficiently rapid to permit the detection of an annual cycle. This probably reflects yearto-year growth rate changes that are not apparent in the U-Th chronology, and that the interval around 4.8 ka BP reflects slower growth than that implied by the relatively low-resolution U-Th chronology. The evolution to gradually more well-developed annual cycles from~4.8 ka BP to~4.2 ka BP (Fig. 12 b-d) potentially reflects steadily increasing growth rates at time-scales not easily discernible using U-Th dating. We therefore suggest that compared with proxies for the ITCZ (the Cariaco Basin Ti% record (light grey) (Haug et al., 2001) with the revised chronology of Kennett et al. (2012)) and the Azores High (the REF-07 stalagmite density record from El Refugio Cave, southern Iberia (red) (Walczak et al., 2015). Also shown is June insolation at 60 N (blue) (Berger and Loutre, 1991). climate leading up to the 4.2 ka event in north Iberia was characterised by gradually decreasing summer rainfall (the increasing bias towards winter rainfall producing lower d 18 O cal values), steadily increasing total annual rainfall (contributing to lower d 18 O and lower Mg values), and a steady increase in GAR-01 growth rates (evident through progressively more clearly expressed annual cycles).
In northern Iberia, we propose that the 4.2 ka BP event was characterised by Mediterranean-like climate conditions that promoted rapid stalagmite growth (e.g., rapid growth is characteristic of stalagmites from Mediterranean-influenced sites like El Refugio cave ( Fig. 1) (Walczak et al., 2015). Similar conditions also existed in the early Holocene, and we suggest that both reflect the summertime influence of the Azores High. The early Holocene pseudo-Mediterranean climate was caused by the insolation-controlled northward migration of the Azores High (Fig. 11), whereas the to 4.265 ka BP (bed) and the YD (e). The model derived % summer rainfall is also shown in (a) (blue line). Note that panel e illustrates a 10-year period during the Younger Dryas interval discussed in Baldini et al. (2015) for comparison with the select mid-Holocene data. The longer YD record is not shown here; see Baldini et al. (2015). Panels cee have identically scaled x-and y-axes and the horizontal dashed line in each panel represents the mean Mg concentration across the entire YD event (11.7e12.880 ka BP) rather than just the 10-year window shown in (e). Vertical grey bars in panel (a) labelled b, c, and d provide the relative timing of the 10-year windows presented in panels b, c, and d, respectively.
climate leading into the 4.2 ka event was caused by the slow migration of the Azores High due to south tracking decreasing insolation. This interpretation is also consistent with observed Northern Hemisphere monsoonal weakening and southward ITCZ migration over the Holocene (Haug et al., 2001) (Fig. 11). The GAR-01 Mg concentration data (Fig. 12) do not contain any evidence for a discrete, short-lived 4.2 ka event; rather the 'event' was the culmination of a~1000-year period of slowly changing climate. The interval including and leading into the 4.2 ka event represents the last prolonged episode of very 'anomalous' climate prior to the emplacement of an essentially modern climate system in northern Iberia.

Conclusions
The results suggest that seasonal rainfall and temperature distribution variability can successfully account for 95% of the record. Only seven timeslices exist where the model does not converge, suggesting the influence of a factor not considered in the model. Interestingly, two of these timeslices are at 13.0 and 12.9 ka BP, coinciding with the YD event's initiation, and may represent poorly understood large-scale atmospheric and oceanic re-organisation associated with this event. The other intervals where the model fails to converge are at 12, 11.4, 9.6, 4.7, and 1.7 ka BP. The modelled seasonality reconstruction can therefore account for the vast majority of the observed d 18 O variability, and presents a~13 ka BP model of seasonal temperature and rainfall shifts in northern Iberia.
The model suggests that the low d 18 O values occurring during the peak YD were primarily the result of low regional air temperatures, consistent with Baldini et al. (2015), but provides new insights into climate seasonality during the event. The model results suggest that the YD was characterised by cold, dry winters and cool, wet summers. Interestingly, this interval of reversed polarity conditions (compared to modern) resulted in d 18 O values that, although low, were higher than if rainfall (or snowfall) had occurred predominantly in the winter. This interpretation is consistent with a variety of research suggesting that the YD was characterised by considerably colder winter but only marginally cooler summer temperatures (Denton et al., 2005); the extreme winter cooling would have also reduced evaporation from the North Atlantic and consequently rainfall during the winter. Extensive winter sea ice across the North Atlantic during very cold YD winters would also have reduced moisture uptake from the sea surface. The presence of equivalent low d 18 O values in the GAR-01 record in both YD and mid-Holocene timeslices is therefore explicable partially by higher than expected values during the YD. Climate in northern Iberia during the early Holocene was affected by a northward shifted Azores High that produced extremely arid summers but sustained winter rainfall. The presence of a Mediterranean-like climate in northern Iberia is consistent with other Iberian and European records tracking both the Azores High and the position of associated westerly winds (e.g., Walczak et al., 2015;Bakke et al., 2009;Wirth et al., 2013). Mid-Holocene climate was similarly dominated by the position and strength of the Azores High as it slowly migrated to the south, tracking insolation. From~7 to~5 ka BP rainfall was strongly biased towards the summer, and the end of this reversed polarity rainfall seasonality interval coincided with the establishment of Mediterranean conditions in southern Iberia at 5.3 ka BP, as well as the shift from a mostly positive NAO to a more negative NAO shortly afterwards (from 5.2 to 4.4 ka BP) (Olsen et al., 2012), again probably due to southward displacement of the Azores High. The model suggests that summers were becoming drier by~4.8 ka BP, a trend culminating in very arid summers at~4.2 ka BP. The lack of summer rainfall and abundant low d 18 O winter rainfall would have biased annually averaged recharge, which is well-integrated at the La Garma Cave site, towards winter rainfall values, consequently explaining the low d 18 O values observed during this time interval.
This peak in summer aridity is probably associated with the globally-observed 4.2 ka event (Booth et al., 2005;Staubwasser et al., 2003). In GAR-01, the event appears as the peak of the summer drying trend that lasted almost 1 kyr, consistent with flowstone data from Italy (Drysdale et al., 2006). This interpretation is supported by a high-resolution LA-ICPMS transect of the 4.2 ka event interval, and we suggest that the 4.2 ka event in northern Iberia was characterised by higher than average annual mean rainfall with very arid summers. A modern polarity in rainfall is achieved at 1.6 ka BP and persists until the present day; the rainfall distribution during the LIA was characterised by a slightly amplified modern rainfall polarity, with drier than modern summers.
We stress however that these modelling results represent one possible scenario, and the geochemical data is potentially explicable by appealing to other mechanisms, such as shifts in moisture mass trajectory, moisture source region, or evapotranspiration regime not detected by the model. Regardless, seasonality shifts are among the most parsimonious explanations for the observed variability during the climatologically subdued Holocene, a conclusion supported by a recent comprehensive review of Iberian climate suggesting that seasonality is a key variable over the Holocene (Moreno et al., 2017). This study therefore presents critical missing information bridging the Holocene and the Pleistocene in a climatologically sensitive region affected by both AMOC and the NAO.
The modelling approach discussed here provides a more detailed view of climate change compared to a semi-quantitative discussion of the d 18 O record in isolation. Although rainfall amount certainly varied from year-to-year across the last 13 ka, the model suggests that seasonality shifts can explain the vast majority of the GAR-01 d 18 O record's variability. This approach potentially provides critical missing seasonality data for northern Iberia and is supported by high resolution trace element data across the 4.2 ka and Younger Dryas events, but other monthly-scale records are needed to confirm the model's output across other intervals. The same modelling approach is also tenable for d 18 O records from other sites and may help highlight aspects of any record that are explicable by appealing to seasonality shifts. Intervals where a model does not converge could highlight climate anomalies forced by unusual factors. Future stalagmite-based palaeoclimate research could benefit from using similar approaches to support interpretations based on geochemical climate proxies alone.