Pattern, style and timing of British–Irish Ice Sheet advance and retreat over the last 45 000 years: evidence from NW Scotland and the adjacent continental shelf

Predicting the future response of ice sheets to climate warming and rising global sea level is important but difficult. This is especially so when fast‐flowing glaciers or ice streams, buffered by ice shelves, are grounded on beds below sea level. What happens when these ice shelves are removed? And how do the ice stream and the surrounding ice sheet respond to the abruptly altered boundary conditions? To address these questions and others we present new geological, geomorphological, geophysical and geochronological data from the ice‐stream‐dominated NW sector of the last British–Irish Ice Sheet (BIIS). The study area covers around 45 000 km2 of NW Scotland and the surrounding continental shelf. Alongside seabed geomorphological mapping and Quaternary sediment analysis, we use a suite of over 100 new absolute ages (including cosmogenic‐nuclide exposure ages, optically stimulated luminescence ages and radiocarbon dates) collected from onshore and offshore, to build a sector‐wide ice‐sheet reconstruction combining all available evidence with Bayesian chronosequence modelling. Using this information we present a detailed assessment of ice‐sheet advance/retreat history, and the glaciological connections between different areas of the NW BIIS sector, at different times during the last glacial cycle. The results show a highly dynamic, partly marine, partly terrestrial, ice‐sheet sector undergoing large size variations in response to sub‐millennial‐scale climatic (Dansgaard–Oeschger) cycles over the last 45 000 years. Superimposed on these trends we identify internally driven instabilities, operating at higher frequency, conditioned by local topographic factors, tidewater dynamics and glaciological feedbacks during deglaciation. Specifically, our new evidence indicates extensive marine‐terminating ice‐sheet glaciation of the NW BIIS sector during Greenland Stadials 12 to 9 – prior to the main ‘Late Weichselian’ ice‐sheet glaciation. After a period of restricted glaciation, in Greenland Interstadials 8 to 6, we find good evidence for rapid renewed ice‐sheet build‐up in NW Scotland, with the Minch ice‐stream terminus reaching the continental shelf edge in Greenland Stadial 5, perhaps only briefly. Deglaciation of the NW sector took place in numerous stages. Several grounding‐zone wedges and moraines on the mid‐ and inner continental shelf attest to significant stabilizations of the ice‐sheet grounding line, or ice margin, during overall retreat in Greenland Stadials 3 and 2, and to the development of ice shelves. NW Lewis was the first substantial present‐day land area to deglaciate, in the first half of Greenland Stadial 3 at a time of globally reduced sea‐level c. 26 kabp, followed by Cape Wrath at c. 24 kabp. The topographic confinement of the Minch straits probably promoted ice‐shelf development in early Greenland Stadial 2, providing the ice stream with additional support and buffering it somewhat from external drivers. However, c. 20–19 kabp, as the grounding‐line migrated into shoreward deepening water, coinciding with a marked change in marine geology and bed strength, the ice stream became unstable. We find that, once underway, grounding‐line retreat proceeded in an uninterrupted fashion with the rapid loss of fronting ice shelves – first in the west, then the east troughs – before eventual glacier stabilization at fjord mouths in NW Scotland by ~17 kabp. Around the same time, ~19–17 kabp, ice‐sheet lobes readvanced into the East Minch – possibly a glaciological response to the marine‐instability‐triggered loss of adjacent ice stream (and/or ice shelf) support in the Minch trough. An independent ice cap on Lewis also experienced margin oscillations during mid‐Greenland Stadial 2, with an ice‐accumulation centre in West Lewis existing into the latter part of Heinrich Stadial 1. Final ice‐sheet deglaciation of NW mainland Scotland was punctuated by at least one other coherent readvance at c. 15.5 kabp, before significant ice‐mass losses thereafter. At the glacial termination, c. 14.5 kabp, glaciers fed outwash sediment to now‐abandoned coastal deltas in NW mainland Scotland around the time of global Meltwater Pulse 1A. Overall, this work on the BIIS NW sector reconstructs a highly dynamic ice‐sheet oscillating in extent and volume for much of the last 45 000 years. Periods of expansive ice‐sheet glaciation dominated by ice‐streaming were interspersed with periods of much more restricted ice‐cap or tidewater/fjordic glaciation. Finally, this work indicates that the role of ice streams in ice‐sheet evolution is complex but mechanistically important throughout the lifetime of an ice sheet – with ice streams contributing to the regulation of ice‐sheet health but also to the acceleration of ice‐sheet demise via marine ice‐sheet instabilities.


Supplementary Information for
Pattern, style and timing of British-Irish Ice Sheet advance and retreat over the last 45,000 years: evidence from NW Scotland and the adjacent continental shelf

Notes and Definitions
Absolute ages -we use the term absolute here to refer to the fact that all ages in this work are on the same timescale and therefore can be compared in an absolute rather than a relative way. We are aware that all Quaternary age assessments are based on models that may be imperfectly constrained.
ka BP and ka cal BP -all ages and dates in the text are expressed in calendar kilo-years before present (ka BP), on a common (or absolute) timescale. In the Results sections, radiocarbon dates are presented in conventional uncalibrated form (in Tables), and subsequently calibrated into calendar kilo-years before present [1950 CE]. For avoidance of doubt, these ages are stated as 'cal ka BP' in the text.
Last glacial cycle -defined here for convenience as Marine Isotope Stages 2-5d (i.e. 116-11.7 ka BP. Broadly equivalent to the Weichselian glacial Stage in chronostratigraphic studies of Northern Europe and offshore regions. We do not include the warm Eemian (N European) = Ipswichian (UK) interglacial stage (MIS 5e) in this definition. Note: MIS = Marine Isotope Stage throughout.
Last Glacial Maximum (LGM) -refers to the global last glacial maximum when eustatic sea levels were at their lowest (100 to 140 m below present) and NH continental ice sheets were generally at their largest (in volume and area). There is no definitive date for this event, however most northern European palaeoglaciological studies now converge on the time period -23-21 ka BP (e.g. Hughes et al., 2016). Some modellers and sea-level scientists prefer a stricter definition centred ca. 24-26 ka BP (Peltier and Fairbanks, 2007) or ca. 22-19 ka BP (Yokoyama et al., 2018) -when eustatic sea level was at its global minimum. A more relaxed time envelope of ~26-19 ka BP or ~30-19 ka BP is preferred by others, especially when synthesising sea level data with northern and southern hemisphere ice-sheet extent information (e.g. Clark et al., 2009;Carlson & Clark, 2012;Yokoyama et al., 2018). The term Local Last Glacial Maximum (LLGM) is sometimes used for specific ice sheets (e.g. Ballantyne & Small, 2018), but this can cause confusion and use of the term 'LLGM' is discouraged.
Greenland Stadials and ice-core chronology -we use the standardised definitions and timings of Greenland Stadials (GS-1, GS-2, etc) and Greenland Interstadials (GI-1, GI-2, etc) throughout this work as recommended by the INTIMATE group (Lowe et al., 2008;Rasmussen et al., 2014). The timescale is defined by the NGRIP isotopic record and uses the Greenland Ice Core Chronology (GICC05), in years before 2000 CE.
Weichselian vs Devensian -we use the term Weichselian throughout this work (referring to glacial Stage), covering the time period from 116.0-11.7 ka BP -starting at the onset of GI-25 and ending at the end of GS-1. Weichselian is now the standard Northern European terminology, used throughout Scandinavia, including Iceland, Faroe and adjacent offshore regions. To avoid confusion and ambiguity we do not use the terms Devensian glacial Stage or Dimlington Stadial, even though these are commonly used by UK-English authors.

Bayesian modelling Methods
Two Bayesian temporal models have been developed in this work: one for the mainland-sourced ice sheet and one for Lewis. Both models were developed using OxCal 4.3 (Bronk Ramsey, 2013). The Bayesian analysis uses Markov chain Monte Carlo (MCMC) sampling to build up a distribution of possible solutions generating a probability, called a posterior density estimate for each sample, these are the product of the prior model and likelihood probabilities. This has produced modelled ages for boundaries (Tables 10 & 11) delimiting several zones or phases. Each phase groups dating information for sites that share a common relationship with all other items in the model. Phases are separated by boundary commands, which generate modelled age-probability distributions for Boundary Limits that bracket events within the retreat sequence (Tables 10 & 11).
The sequence model was run in an outlier mode to assess outliers in time using a student's tdistribution (p < 0.05) to describe the distribution of outliers with a scaling of 10-10,000 years (Bronk Ramsey, 2009). Complete outliers were given a probability scaling of p = 1, with ages showing a poor model fit given values of p < 0.25, p < 0.5, p < 0.75 and p < 0.90 on a scale of increasing severity. Outliers here, in part, reflect the site-specific treatments (see Results) but also assess the fit of individual age measurements in the different phases of grouped dating information, as well as their fit in the overall retreat sequence. Chronological measurements showing a strong fit to the model have individual agreement indices >60 (Bronk Ramsey, 2009, 2013 Fig. 32). Ultimately, Bayesian analysis produced a conformable age model with an overall agreement index of 158% for the Minch Ice Stream and 138% for Lewis, exceeding the 60% threshold advocated by Bronk Ramsey (2009,2013). Note that ages relating to ice-sheet fluctuations before ~35 ka BP were not included in the Bayesian modelling for reasons of computational complexity. Figure S1: Sedimentology, geochronology and geophysical properties of seabed core JC123-035VC. From left to right: X-radiograph, interpreted lithofacies boundaries, lithological log, main lithofacies codes, P-wave velocity, gamma-ray attenuation (density), magnetic susceptibility and undrained shear strength values all plotted on same depth scale. Ages of AMS-dated samples given in radiocarbon years (italic) and calibrated yrs BP (bold font) (see Table 8). Geophysical properties data measured using a Geotek MSCL-S at 2 cm intervals. Gaps are missing data (end caps). Figure S2: Sedimentology, geochronology and geophysical properties of seabed core JC123-036VC. From left to right: X-radiograph, interpreted lithofacies boundaries, lithological log, main lithofacies codes, P-wave velocity, electrical resistivity, gamma-ray attenuation (density) and magnetic susceptibility values all plotted on same depth scale. Ages of AMS-dated samples given in radiocarbon years (italic) and calibrated yrs BP (bold font) (see Table 8). Geophysical properties data measured using a Geotek MSCL-S at 2 cm intervals. Gaps are missing data (end caps). Figure S3: Sedimentology, geochronology and geophysical properties of seabed core JC123-044VC. From left to right: positive X-radiograph, interpreted lithofacies boundaries, lithological log, main lithofacies codes, P-wave velocity, electrical resistivity, gamma-ray attenuation (density) and magnetic susceptibility values all plotted on same depth scale. Ages of AMS-dated samples given in radiocarbon years (italic) and calibrated years BP (bold font) with uncertainties (see Table 8). Geophysical properties data measured using a Geotek MSCL-S at 2 cm intervals. Gaps are missing data (end caps). Note: see text for explanation and relevance of chronology Figure S4: Sedimentology, geochronology and geophysical properties of seabed cores JC123-018VC. From left to right: positive X-radiograph, interpreted lithofacies boundaries, lithological log, main lithofacies codes, P-wave velocity, gamma-ray attenuation (density) and magnetic susceptibility values all plotted on same depth scale. Age of AMS-dated sample in radiocarbon years (italic) and calibrated years BP (bold font) with uncertainties (see Table 8). Geophysical properties data measured using a Geotek MSCL-S at 2 cm intervals. Gaps are missing data (end caps).  positive X-radiograph, interpreted lithofacies boundaries, lithological log, main lithofacies codes, Pwave velocity, gamma-ray attenuation (density) and magnetic susceptibility values all plotted on same depth scale. Geophysical properties data measured using a Geotek MSCL-S at 2 cm intervals. Gaps are missing data (end caps). Figure S7: (above) sedimentology, geophysical properties and acoustic stratigraphy of seabed core JC123-029PC, Sound of Raasay, Inner Minch. From left to right: positive X-radiograph, interpreted lithofacies boundaries, lithological log, main lithofacies codes, magnetic susceptibility and P-wave velocity values all plotted on same depth scale. Ages of AMS-dated sample in radiocarbon years (italic) and calibrated years BP (bold font) with uncertainties (see Table 8). Geophysical properties data measured using a Geotek MSCL-S at 2 cm intervals. Gaps are missing data (end caps). Note the break between 80-300 cm to allow for Figure to be scaled appropriately on printed page. Table S1.