Exploring the extent to which fluctuations in ice‐rafted debris reflect mass changes in the source ice sheet: a model–observation comparison using the last British–Irish Ice Sheet

The British and Irish Ice Sheet (BIIS) was highly dynamic during the Late Quaternary, with considerable regional differences in the timing and extent of its change. This was reflected in equally variable offshore ice‐rafted debris (IRD) records. Here we reconcile these two records using the FRUGAL intermediate complexity iceberg–climate model, with varying BIIS catchment‐level iceberg fluxes, to simulate change in IRD origin and magnitude along the western European margin at 1000‐year time steps during the height of the last BIIS glaciation (31–6 ka bp). This modelled IRD variability is compared with existing IRD records from the deep ocean at five cores along this margin. There is general agreement of the temporal and spatial IRD variability between observations and model through this period. The Porcupine Bank off northwestern Ireland was confirmed by the modelling as a major dividing line between sites possessing exclusively northern or southern source regions for offshore IRD. During Heinrich events 1 and 2, the cores show evidence of a proportion of North American IRD, more particularly to the south of the British Isles. Modelling supports this southern bias for likely Heinrich impact, but also suggests North American IRD will only reach the British margin in unusual circumstances.


Introduction
An ice sheet with a large marine component, such as the West Antarctic Ice Sheet today, is suspected to be dynamic, with major change occurring quickly but not spatially synchronously, leading to major consequences for iceberg flux rates, ocean circulation and, potentially, sea level (Hughes, 1975;MacAyeal, 1992;Bamber et al., 2009;Thomas, 2014). The best studied past analogue for such a dynamic marine glacial system is the British and Irish Ice Sheet (BIIS) of the last glacial period (Clark et al., 2018). The evolution, and particularly retreat pattern, of the last glacial BIIS has been studied from both a largely 'deep ocean' perspective, driven by temporal and spatial distributions in along-margin ice-rafted debris (IRD) flux, sea surface temperature (SST) and lithology (e.g. Scourse et al., 2009); and a largely 'terrestrial and shelf' perspective, driven by the spatial and temporal distribution of glaciogenic landforms (e.g. Clark et al., 2012Clark et al., , 2020. Both of these approaches lead to the conclusion that the BIIS was indeed dynamic in terms of the asynchronous nature of the peak advance and retreat patterns. However, what is currently lacking in the literature for any past ice sheet, is a formal reconciliation of the detailed reconstruction of these two perspectives. Our approach to this reconciliation is to use the wellconstrained BRITICE-CHRONO reconstruction of the growth and retreat of the British-Irish Ice Sheet during the last glacial (31 to 15 ka BP) to compute potential iceberg discharges into the ocean (Clark et al., 2020). To track the transport and deposition of the resulting IRD we use an intermediate complexity coupled iceberg-climate model (Levine and Bigg, 2008) to make IRD predictions across this time period, accounting for temporal climate and ice sheet change. These model ice sheet-iceberg-ocean predictions are then compared with the empirical record of IRD at five well-described locations. Put differently, the novelty here is that we use iceberg trajectory modelling to permit exploration and reconciliation of the regional ice sheet and IRD records. To our knowledge this is the first time such an approach has been undertaken; the significance being that a large body of literature (e.g. Hemming, 2004;Peck et al., 2007;Hall et al., 2011) implicitly assumes that there is a tight, predictable and easily interpreted relationship between ice mass change and IRD occurrence. We aim to assess how simple or complex this relationship is using a well-constrained ice sheet with a rich IRD record.
The BIIS context for this analysis is first described, followed by a discussion of the various IRD and iceberg flux data, as well as the climate model to be used. A simulation of the evolution of the combined flux and iceberg distribution, using 1000-year time slices from 31 to 16 ka BP is then presented (the ice sheet had disappeared by 15 ka BP). The sensitivity of the simulation to both Heinrich event intervention from the western Atlantic, and two different reconstructions for the ice flux from the Porcupine Bank, off western Ireland, is then tested. The paper concludes with a discussion of the results and consideration of the consequences for the terrestrially focused reconstruction.

BIIS context
We focus our investigation within the BIIS context around a series of questions to gain insights into the interactions: 1. Does IRD typically come from the obvious local source, as suggested by Scourse et al. (2009), or do ocean currents confuse the picture? 2. To what extent can icebergs originating from other ice sheets affect regional records? 3. Do IRD events in the marine record reflect mean iceberg fluxes, or mostly record advance or retreat phases? 4. To what extent are the core sites in, rather than avoiding, iceberg alleys, and is this a significant factor that should be acknowledged in making IRD interpretations? 5. Do variations in ocean currents and the position of the North Atlantic Polar Front (NAPF) over time either directly or indirectly affect the transport of icebergs?
The investigation is driven by the west European offshore IRD perspective epitomised by Scourse et al. (2009) and the comprehensive BIIS palaeo-ice sheet reconstruction from the BRITICE-CHRONO project (http://www.britice-chrono.group. shef.ac.uk/; reported in Clark et al., 2020). The reconstructed palaeo-ice stream routes from BRITICE-CHRONO are shown in Fig. 1, along with the IRD marine core positions. Note that ice streams internal to the continental shelf, as well as from the independent Shetland ice mass, are not considered in our investigation.
The open ocean IRD perspective of the late BIIS evolution appears to be simplified by the vast majority of the offshore IRD originating from the BIIS itself, where IRD lithologies can be fingerprinted (Peck et al., 2007;Scourse et al., 2009). The main possible contaminant of IRD from the Laurentide Ice Sheet (LIS) of North America is only found in the marine record along a small part of the central west European margin, off southwest Ireland, for short intervals, during Heinrich events (Fig. 2;Scourse et al., 2009). It is possible that IRD sourced from Greenland, Iceland or the Faroes has been conflated with BIIS-sourced IRD in this analysis. However, the observed lack of LIS material off the glacial British Isles is consistent with the possible trajectories from the iceberg modelling work of Bigg et al. (2010). The first two questions above are driven by these findings from Scourse et al. (2009). One major question raised by their IRD analysis was the indicative meaning of IRD: does an increase in IRD flux represent an extension in the calving margin during ice sheet growth, or rapid collapse during deglaciation, or both (Question 3 above)?
Questions 4 and 5 above are driven by how one expects the marine IRD record to be driven by changes in climate, and associated changes in the geography of ice sheet extent and thickness. One argument is that the BIIS will expand or reorientate southwards during colder stadials, resulting in IRD expanding in the southern shelf region as more ice is calved to the south, while the IRD record north of the Polar Front will experience a period of constant IRD flux during these southward excursions (Scourse et al., 2009). During our period of interest there are Greenland interstadials at~23 (GIS2), 27 (GIS3) and 28.5 (GIS4) ka BP (Rasmussen et al., 2014). The combination of this hypothesised mechanism and the IRD record of Fig. 2 led Scourse et al. (2009) to suggest early, major advance of shelf ice off western Scotland after GIS4, being maintained until around GIS2. A major southward expansion of ice is consistent with the stadial between GIS2 and GIS3, with the central shelf around Porcupine Bank maintaining significant IRD flux throughout, being near or north of the NAPF throughout much of the period.
The terrestrial and shelf-centred view of the evolution of the peak and retreat of the BIIS, based on Clark et al. (2012), updated and consolidated by Clark et al. (2020) in a much-enhanced datarich and ice sheet modelling analysis, starts from evidence that the BIIS had a significant ice core within northern Britain and Ireland  Scourse et al. (2009), and the main ice streams from the BRITICE-CHRONO reconstruction of Clark et al. (2020). The topography is taken from Weatherall et al. (2015). [Color figure can be viewed at wileyonlinelibrary.com] around 30 ka BP, but not extending across the North Sea or into southern Ireland. A rapid and major ice sheet expansion is thought to have occurred around 27 ka BP, with ice extending out to the continental shelf edge around much of western Britain and Ireland (Peters et al., 2016;Callard et al., 2018). This is reconstructed to be short-lived, with significant retreat from the western Scottish shelf, around the time of, or just after, GIS3 (Ó Cofaigh et al., 2019). This was followed by another period of rapid expansion (25-24 ka BP), although at this time centred in the Irish Sea and extending towards the Isles of Scilly (Smedley et al., 2017;Scourse et al., 2019). This advance is reconstructed to be short-lived also, but the ice extent, and calving fluxes, from the core BIIS remained similar for several thousand years. It was then not until after 20 ka BP that the main BIIS began to retreat back onshore, shrinking gradually back to its Scottish core by 16 ka BP.
The reasons for both of these reconstructions given above are partially ascribed to climate change; for example, through high relative sea level (Ó Cofaigh et al., 2019), but also to internal ice stream dynamics (e.g. Sejrup et al., 2016;Gandy et al., 2018). Here, through use of an intermediate complexity climate and iceberg model (Levine and Bigg, 2008), the linkage between the two perspectives is studied by attempting to reconcile the temporal and spatial variability of both the iceberg calving flux and the offshore IRD flux. As part of this reconciliation the question about the sensitivity of the marine core locations will also be addressed; this is Question 4 as posed earlier.

IRD data
The main IRD data which this paper uses are described in Scourse et al. (2009) and Table 1. There are five cores spanning the west European shelf from northwest Scotland to Brittany. The core locations are shown in Fig. 1, and the IRD data in Fig. 2. As shown in Table 1, the various cores cover different age limits, but here interest is restricted to 31-16 ka BP. Note that we only use core MD04-2829CQ (Hall et al., 2011) from Rosemary Bank and not the other in this location discussed in Scourse et al. 2009, namely DAPC-2, as the latter does not cover the full extent of the time period under consideration.
As can be seen from Table 1, the cores have varying sedimentation rates so age models for all cores are based on tuning the percentage frequency of the >150 μm fraction of polar planktonic foraminfer Neogloboquadrina pachyderma sinistral (Nps) to the GISP2 ice core δ 18 O record (Scourse et al., 2009). The rationale for this tuning is that the %Nps is a proxy for the meridional movement of the NAPF, and so linked to northern Atlantic climate variability, as will be the Greenland ice sheet's properties. A number of 14 C and tephra ages help constrain these records (for more details, see Scourse et al. (2009) and papers referenced therein). The one exception to this procedure is core MD95-2006, from the Barra Fan. While prior to c. 30 ka BP the sedimentation at this core seems largely hemipelagic, and so a  Scourse et al. (2009). Note, however, that Barra Fan shows the concentration, rather than flux, with units >250 μm grains g −1 (see text for discussion). Cores shown from north to south (when reading left to right, top to bottom; see Table 1 for core names and positions). Note that the timing of H1~15-16 ka BP, H2~24-25 ka BP and H3 is~30-31 ka BP. tuning approach can be used (but for the >250 μm fraction), after c. 30 ka BP the location of this core on an active trough mouth fan led to considerable turbiditic sedimentation in its record (Knutz et al., 2001). This makes interpretation of the details of the IRD record more difficult. There are some 14 C dates in this section to help constrain the age model (Knutz et al., 2001); these agree with the reconstructions presented in the next section. The Barra Fan IRD record is therefore discussed here, but requires care in its interpretation.
For interpretation of the IRD record it is important to note that provenance analysis carried out in Scourse et al. (2009), and papers referenced therein for specific cores, showed that most of the IRD record comes from locally sourced lithic material (see Peck et al. (2007) for more details). This is likely to have been largely carried by icebergs derived from the catchment areas of adjacent ice streams. While it is very difficult to discriminate IRD derived from a marginal sea-ice source as opposed to that transported by icebergs (St. John et al., 2015), the former is unlikely to be widespread along the eastern Atlantic margin and will not be transported far before melting (Vettoretti and Peltier, 2013). The only long-distance transport of non-BIIS IRD found in Scourse et al. (2009) was limited material carried from the LIS during Heinrich events, specifically H1 and H2 in our time period. This is only present in significant quantities at those times in the Goban Spur core, with a marginal H2 signal in the Porcupine Bank core. Note that the absence of LIS dolomitic carbon around H3, at 30 ka BP, is consistent with this event originating from a different, possibly Arctic, ice sheet (Bigg et al., 2011).

Ice volume calculations
To estimate the ice flux from the BIIS between 31 and 16 ka BP, steady-state ice sheet simulations forced to coincide with the reconstructed isochrones of Clark et al. (2020) were performed at 1000-year intervals. Ice thickness was derived using an ice sheet model which assumes perfectly plastic ice flow (Gowan et al., 2016). This plastic ice sheet model requires three inputs: ice extent, basal shear stress and basal topography. Ice extent was prescribed to the plastic ice sheet model from the reconstruction of Clark et al. (2020). Basal shear stress was assigned according to glacial geomorphological units using the map of Gandy et al. (2018). The general bathymetric chart of the oceans (www.gebco.net; Weatherall et al., 2015) was used as an input basal topography. All simulations had a horizontal resolution of 5 km.
A balance flux approach, whereby the flux of ice at a point is assumed to balance upstream accumulation, was used to estimate ice flux through the major outlets of the BIIS (Le Brocq et al., 2006). The balance flux approach has been shown to reasonably recreate the ice flux patterns of the Antarctic ice sheet (Le Brocq et al., 2006). The balance flux algorithm requires inputs of ice thickness (derived from the plastic ice sheet model) and accumulation. Unfortunately, we do not have an accumulation proxy or climate-ice sheet model which is consistent with the reconstruction of Clark et al. (2020). Therefore, variation in latitudinal and longitudinal gradients of accumulation were assumed to be the same as for contemporary precipitation, an approach which has been shown to be successful for modelling aspects of the BIIS (Patton et al., 2016). Ice flux was then integrated from the balance flux calculations across the margin of each of the major BIIS outlets and the Norwegian Channel Ice Stream (Fig. 1). The volume of ice in each catchment was calculated for every time interval. Changes in volume (either a growth during ice advance, or reduction during ice retreat) were distributed to adjacent time steps to correct for ice volume storage or loss in the final estimates of iceberg flux, as shown in Fig. 3. Though both the plastic model and balance flux algorithm contain large simplifications of the mechanics of ice sheet flow and palaeoclimatic conditions, the derived iceberg fluxes are of a similar order of magnitude to previous estimates  and provide a first-order estimate of the importance of different catchments for iceberg flux through  (Fig. 3). These simple models match the reconstructed ice limits and deglaciation chronology of the BIIS in a manner that more complex models may struggle to achieve.

FRUGAL intermediate complexity model
The FRUGAL (Fine ResolUtion Greenland And Labrador) intermediate complexity global climate model includes coupling between ocean, radiative-advective atmospheric, simple advective-thermodynamic sea-ice and iceberg trajectory models (Levine and Bigg, 2008). FRUGAL has been used in a number of palaeoclimate studies over a range of times including the last glacial period (Levine and Bigg, 2008;Bigg et al., 2010Bigg et al., , 2011, the penultimate glaciation (Green et al., 2010) and the Early Pleistocene (Rea et al., 2018). The model uses a curvilinear grid with the North Pole in Greenland, giving an enhanced resolution in the Arctic and North Atlantic (Wadley and Bigg, 2000). Here we use a fine resolution grid of 182 ×211 cells, which is equivalent to approximately 2°longitude by 1.5°latitude in the Southern Hemisphere, but with a resolution of 20 km around the Greenland coast. The details of the ocean and sea-ice model configuration are given in the similar resolution, but present-day, study by Wilton et al. (2015). The atmospheric part of FRUGAL is a simple radiative-advective atmosphere (Fanning and Weaver, 1996) that allows for advection of water vapour. All the simulations discussed in this paper use a monthly varying glacial wind stress, which has no feedback from the SST field (Levine and Bigg, 2008). The basic topography and bathymetry of the model is as described in the last glacial FRUGAL simulations (e.g. Levine and Bigg, 2008), but with minor modifications around the British Isles to reflect knowledge of ice sheet extent (Clark et al., 2012), and taking into account the finer resolution of the model grid in the shaping of local bathymetry for important gateways and shallow seas. The iceberg module has both dynamic and thermodynamic components, although in this study there is no feedback between the iceberg model and the rest of the FRUGAL system. The dynamical processes included in the iceberg model are ocean/atmosphere/ sea-ice drags, the Coriolis force, pressure gradients in the surrounding ocean, and wave radiation. The main thermodynamic processes represented are basal melting, buoyant convection, wave erosion, and several smaller terms due to sublimation and latent heat transfer. Icebergs may roll over (Wagner et al., 2017) and grounded icebergs are allowed to melt instantaneously (Levine and Bigg, 2008). Model icebergs are divided into 10 different size classes, ranging from 0.491 to 492 × 10 9 kg in mass, based on observations of present-day Arctic and Southern Ocean icebergs, excluding giant icebergs. Each model berg is assigned a scale factor appropriate to that size of berg from its specific seed site. Summing the mass of each berg, multiplied by its scale factor, gives the average ice discharge per quarter year expected from that seed site (see Levine and Bigg (2008) for a depiction of this in the last glacial period). A 1000-year spin-up simulation was run, using orbital parameter and atmospheric CO 2 levels from 21 ka BP (as in Levine and Bigg (2008)). By the end of this simulation all largescale ocean and atmospheric fields had reached an essential state of equilibrium. A set of sensitivity experiments were then carried out. In this set, 16 × 10 year simulations were carried out, with each simulation using a different orbital parameter configuration with steps of 1000 years in the parameters, as determined from the equations in Berger (1988), from 31 ka BP to 16 ka BP. The climate changes for each step, being due to minor changes in eccentricity, obliquity and precession (Fig.  S1), are small and the equilibrium state has shifted smoothly within a few years (Fig. S2). Tests where atmospheric CO 2 was changed similarly with time, in line with ice core changes (Lüthi et al., 2008), showed minimal difference across the time period studied. In the set of 10-year simulations the volume change iceberg flux (Fig. 3) is seeded around the British Isles, with iceberg fluxes elsewhere globally, including around the North Atlantic, as given in Levine and Bigg (2008).
A sensitivity study to investigate the local impact of adjustment to iceberg release sites to take account of different scenarios for ice extent around the Porcupine Bank (Peters et al., 2016) was also carried out. To investigate fully the impact of Heinrich events on the possible supply of LIS IRD to the British margin, 10-year simulations starting from year 5900 of the St. Lawrence and Hudson Strait Heinrich event studies of Levine and Bigg (2008) were run, but including giant icebergs of 20 and 40 km in length in the iceberg seed file, replacing the two largest iceberg classes of the original scheme, but preserving total iceberg flux. The intention of this test was to address Question 2 by ascertaining whether icebergs more typical of those released from LIS ice sheet collapse would reach further east before melting.

Ocean properties
While there are minor changes in the global ocean and atmospheric properties between each simulation across the set of 16 time slices, from 31 to 16 ka BP (Table S1), the key features of North Atlantic Ocean currents and SST remain the same. The typical 30 m ocean current is shown in Fig. 4, with an expanded view for the northeast Atlantic off western Europe shown for a selection of time slices in Fig. S3. The basic ocean circulation consists of a sub-tropical gyre, separating from the  North American coast off New England, and sub-polar gyres in the northwestern Atlantic and Norwegian-Greenland Sea. Note that the five IRD sites, whose positions are marked on Fig. 4, are well placed to capture IRD in the coastal currents off glacial Britain, and to show north-south movement of the NAPF boundary. This boundary, separating the southward return flow of the sub-tropical gyre from the northward flow along the British Isles entering the sub-polar gyre systems, occurs north of Ireland, near the Donegal and Malin-Hebrides ice streams (Fig. 1). This is the case for all time slices (Fig. S3), although the similarity is a reflection of the wind stress being the same in all time slices, even if temperatures and salinities differ because of the orbital parameter variation. The Greenland-Norwegian Sea sub-polar gyre is driven by a European coastal shelf current, with the western outflow being both through the Denmark Strait, but also southeast of Iceland. The northwestern Atlantic sub-polar gyre is centred in the Labrador Sea, but there is a general weak, northwestward circulation in much of the northern Atlantic (Figs. 4 and S3). Note that this tends to be weakest offshore from the BIIS, a feature that we will return to in the discussion.

Iceberg distributions
The simulated iceberg density field across the northern Atlantic is shown in Fig. 5 for a selection of time slices, to illustrate the impact of the temporal evolution of BIIS iceberg fluxes (Fig. 3). Note that the seeded iceberg flux from all other ice sources outside the BIIS is fixed for all time slices, and what is shown is the average density for the last two years of a 10-year iceberg run, to ensure the simulation's near-surface properties have reached equilibrium with each change in orbital parameters across the set. The variation in Fig. 5 is therefore from a mix of forced BIIS iceberg release changes, changes in iceberg trajectories and lifetimes due to simulation-dependent ocean temperature and currents, and random noise due to the two-year data window.
At 31 ka BP, the earliest time slice shown in Fig. 5, the only source of BIIS icebergs is from the Minch Ice Stream (Fig. 3), so there are very few icebergs present off the British Isles. The simulations overall show that icebergs from other, non-BIIS, sources are rare and only appear at a few time slices, and exclusively from Scandinavia. At the IRD sites, undetectable levels of IRD are found for this time slice (Fig. 6), consistent with Fig. 5. The low iceberg density area off northwestern Iberia comes from North American icebergs; this is clearer in later time slices.
By 27 ka BP the BIIS had expanded significantly in all directions, with marine ice spreading out to the continental shelf edge around much of western Britain and Ireland (Fig. 3). This led to icebergs being calved into the ocean extensively, but in the model they tend to be trapped in the near-shelf area (Fig. 5), with transport southwards towards Iberia and northwards into the Norwegian Sea, but few spreading westward into the Atlantic proper. The model suggests there is a strong source separation for the icebergs that feed the IRD sources, with those cores north of the NAPF being supplied largely from the Malin-Hebrides ice stream, while those south of the NAPF are fed from the Porcupine Bank, and further south (the Armorican Fan), from the Irish Sea (Fig. 6).
Over the next few thousand years, while the iceberg flux stabilises and calving retreats from the continental shelf edge in places, the basic modelled iceberg distribution, both in spread (Fig. 5) and direction (Fig. 6), remains similar. It is only after 21 ka BP that the retreat of the BIIS from its marine margins intensifies, so that by the end of our study period, 16 ka BP, only iceberg calving off Moray Firth and the Malin-Hebrides ice stream continues in the north, with a very limited Irish Sea flux, leaving only a patchy iceberg distribution around the British Isles (Fig. 5). Note that the source regions of IRD remain essentially the same, with only the northernmost two locations retaining modelled IRD signatures by 16 ka BP (from Malin-Hebrides; Fig. 6).

Sensitivity tests
Two sensitivity tests were carried out to test hypotheses about alternative iceberg sources. Firstly, there is some argument over the extent of the ice bridge out to Porcupine Bank during the peak glacial period (Peters et al., 2016). Did the ice extend all the way to Porcupine Bank, essentially leading to a northern and southern calving source from the resulting ice bridge, or was there a less extensive, single source ice stream in the area? Simulations were therefore re-run, with the Galway Porcupine iceberg flux split equally between a northern and southern seed site (Fig. 7). To better simulate the rapid retreat of the Irish Sea Ice Stream the Irish Sea seed point was also moved towards the coast (except at 26 ka BP; Fig. 7). As with the previous runs, no Galway Porcupine icebergs moved north, and the overwhelmingly dominant source of icebergs from Galway Porcupine south was from the southern release site. Only a very few northern Galway Porcupine source icebergs contributed to the Porcupine Bank and Armorican Fan IRD sites (Fig. S4).
The other sensitivity run involved testing the hypothesis that production of giant icebergs from a Heinrich event collapse of the LIS would enable some LIS bergs to survive the trans-Atlantic crossing. We have already seen that modelling ordinary icebergs leaving the LIS means they do not reach the British Isles (Fig. 4). Using a Heinrich event simulation, where ocean conditions are modified by the much larger LIS iceberg flux of a Heinrich event (Levine and Bigg, 2008), and running simulations with both ordinary and giant icebergs forming part of the flotilla entering the Atlantic, means that a few icebergs, in both cases, are modelled to reach southern Britain (Fig. 8).

Comparison of model and IRD data
The temporal resolution of the IRD record in the five cores is variable, and the model simulations are only for each 1000 years Figure 6. Bar charts of the ice-rafted debris (IRD) source variation over time within 2°of the 5 IRD marine cores, using modelled iceberg mass (in Tg) as a proxy (see Table 1 for core names and locations). Note that the dotted line is the total BIIS iceberg flux (in km 3 yr -1 ; see right-hand scale), as reconstructed in Fig. 3. See Figs. 1 and 3 for release sites and fluxes, respectively. [Color figure can be viewed at wileyonlinelibrary.com] between 31 and 16 ka BP. The units of the two iceberg measures are also different. The cores record IRD in grains cm -2 kyr -1 , in a size range (grains g -1 in the case of Barra Fan; see Fig. 2), while the model results give the mass of icebergs found within 2°of a core site over a two-year period (Fig. 6). For each core, Fig. 9 therefore shows both records averaged over 1000 years across 31-16 ka BP, with the respective units on opposite y-axes. It is the correlation of the two sets of data, and particularly the timing of major respective peaks, that is of interest. Both IRD sourcing by Scourse et al. (2009) and the identification of iceberg origin by the model suggest local sourcing of the majority of material, with northward movement north of the Porcupine Bank and southward movement from here and south of this area throughout the whole period. Thus Fig. 9 shows that, for the BIIS, Question 1 is answered in the affirmative, with local sources dominating the IRD deposition. This is also found to be the case in other areas close to past ice sheets, such as the Denmark Strait (Andrews et al., 2014) and the Labrador Sea (Pearce et al., 2015). However, IRD is well known to be able to travel a very long way in the right oceanographic situations (e.g. Becquey and Gersonde, 2002;Hemming, 2004;Bigg et al., 2010).
A key question for this paper (Question 3) is to what extent the IRD fluxes represent balance ice flux, advance or retreat phases of the ice sheet providing the IRD. The iceberg-climate model was forced by assumed balance ice fluxes. At each core site the Fig. 9 comparison shows that the model captures the peak IRD flux around 25-27 ka BP, although in some places with a 1000-year mismatch, perhaps reflecting the resolution of the ice-sheet reconstruction. However, the details of the comparison vary. These will be examined in turn, from north to south, starting with Rosemary Bank (see Fig. 2 for the detailed IRD record and Fig. 9 for the smoothed comparison). Here the IRD observations show the presence of IRD throughout the period of interest, with a sharp peak around 28 ka BP, followed by a more prolonged peak~24-25.5 ka BP. The iceberg model also shows significant background levels of IRD throughout, but with one peak, around 27 ka BP. This peak occurs because of the strong peak in local catchment balance at this time (Clark et al., 2020). Past reconstruction has suggested a rapid retreat from the continental edge, but then stabilisation (Clark et al., 2012(Clark et al., , 2020 inshore. It is possible that the later (24-25 ka BP) observed IRD peak is associated with the collapse of the shelf-edge extension of the BIIS, rather than reflecting a period of peak expansion of the ice sheet. This has been hypothesised for other areas (Andreassen et al., 2014) and is a classic interpretation of Heinrich event IRD peaks (Hemming, 2004). Such rapid retreat would be smoothed out over the 1000-year time steps of the ice stream flux model and so not be seen as an abrupt step in the iceberg model forcing.
Barra Fan is the next site south. This core only has an IRD record until c. 22 ka BP because of the strongly turbiditic nature of the sediments. However, prior to this time there is excellent agreement in the timing of the peaks in both IRD flux and modelled iceberg number, and some similarity in magnitudes, with the earlier peak being reduced relative to the two later peaks. The fairly continuous ice-streaming activity suggested by Callard et al. (2018) is borne out by both model and IRD observations. The mid-margin core at Porcupine Bank showed one main peak in both IRD and iceberg model data. The iceberg peak is~1000 years later than the IRD peak, but this difference may be due to the details of the core age model or the ice sheet-wide reconstruction. Again, both suggest a period of flux build-up for a few thousand years prior to the peak, from an initially low level. The modelled iceberg number (Fig. 9) is maintained at a higher level for longer after the peak than is the case for the IRD signal. The rapid decline and then re-establishment of IRD levels after the 24-25 ka BP peak (Fig. 2) suggests that the actual ice sheet can respond very rapidly (Peters et al., 2016); faster than the time stepping of either the model or ice flux calculations permit.
The penultimate core in this sequence is Goban Spur, off southwestern Ireland. Both model and IRD observations suggest a persistent low background level of iceberg presence here, with both possessing one major peak. The actual (IRD) peak occurs 1-2000 years after that in the iceberg model. Given the sharp timing of the Irish Sea extension in the balance flux calculations (Fig. 3), yet the lack of an early IRD peak either at the Goban Spur, or on the Armorican Fan around 27 ka BP, suggests that the reconstructions of peak extent around 25 ka BP (Ballantyne et al., 2017;Scourse et al., 2019) are more consistent with the IRD record. The final core IRD record, at the Armorican Fan, has both notable similarities and differences with the modelled iceberg fluxes. The lack of an observed 27 ka BP peak has already been mentioned. However, both records capture a peak around 24 ka BP, as the Celtic Sea advance rapidly retreated (Scourse et al., 2019). Both show some IRD more recently in the record, but the modelled fluxes are much greater. A likely reason for this difference is because the ocean and iceberg model coastline does not alter during the return of the Celtic and Irish Seas (Fig. 7), meaning the model releases icebergs too far out to sea after~25 ka BP. Model icebergs thus reach the Armorican Fan from the Irish Sea (see Fig. 6) when in reality they either had much further to travel, thus melting, or were grounded within the contemporary shallow water of the Irish and Celtic Seas.
From consideration of each of the IRD records, and comparison with the iceberg-climate model's iceberg density, it is clear that the answer to Question 3 is that, while prolonged ice stability and so slowly evolving balance ice fluxes lead to a sustained marine IRD record, peaks in this record relate strongly to change, being particularly sensitive to rapid retreat. Nevertheless, rapid advance can also lead to IRD peaks in circumstances where the advance allows icebergs to enter the open ocean.

The influence of the LIS in the British margin IRD record
It is clear from the provenance analysis of the IRD in the five cores examined here that the only periods of notable presence of sediment from the LIS are during Heinrich events 1 and 2 (Scourse et al., 2009). Even during those major climatic events, LIS material is only found on the Goban Spur and at very small levels on the Porcupine Bank; the latter mostly during H2. This agrees with the iceberg model results. During 'normal' periods, as shown in Fig. 5, modelled icebergs from the LIS may just Figure 9. Comparison of the observed ice-rafted debris flux (blue; in >150 μm grains cm −2 kyr -1 )/1000, except for Barra Fan, where units are >250 μm grains g −1 ) with the modelled iceberg mass flux (orange; Tg within 2°of the core site, over years 9-10 of the simulation) at the five core sites (see Fig. 1). Note, both records have been smoothed with a 1000-year averaging, for comparability. Cores shown from north to south (when reading left to right, top to bottom; see Table 1 for core names and positions). [Color figure can be viewed at wileyonlinelibrary.com] reach the Iberian margin at low concentrations, but they do not penetrate further north. The prevailing surface currents act against any northward movement along the European margin north of Iberia (Fig. 4), preventing LIS icebergs reaching the British margin. However, during Heinrich events, the cooling of the climate and ocean, combined with the impact of a colder, fresher North Atlantic on the ocean currents (Levine and Bigg, 2008), means that some icebergs can penetrate to the European margin off southwestern Ireland (Fig. 8). These do not have to be giant icebergs from an ice shelf collapse, but the presence of such icebergs during a Heinrich event (Hemming, 2004), combined with the sustained period of iceberg supply, will enhance the IRD signal from LIS sources at the Goban Spur during Heinrich events. The presence of limited LIS IRD at Porcupine Bank during H2 suggests the Polar Front was a little further south then, and some icebergs from the Hudson Strait managed to remain north of this boundary in crossing the Atlantic (Hemming, 2004). This more southerly position of the Polar Front during Heinrich events is consistent with the Nps data from this and other sites (Scourse et al., 2009;Haapaniemi et al., 2010). Both the IRD and model record are therefore consistent, with the answer to Question 2 being that only exceptional distal ice sheet change affects IRD levels off the BIIS.

Reconciling model and core ocean properties
It has just been seen that the Nps levels (indicative of cold polar waters) off the British margin during Heinrich events are consistent with previous model results for such conditions, where the East Atlantic is colder (Levine and Bigg, 2008;Bigg et al., 2011) and IRD from the LIS can reach the southern section of our study area. Nevertheless, there are some inconsistencies in the background Nps signal found by Scourse et al. (2009) along the western European margin. As expected, background levels at Porcupine Bank (50-60%) are lower than those at the more northerly Rosemary Bank (60-80%). Levels at the southernmost site of MD95-2002 are even lower during the first half of our study period (10-40%). These all fit with the model's Polar Front typically being off western Ireland. However, %Nps levels at Goban Spur, off the Isles of Scilly (Fig. 1), are 80-90% for much of our study period; MD95-2002 %Nps levels approach 100% during c. 19-15 ka BP (Scourse et al., 2009). Neither of the more northern cores displayed dramatic increases in %Nps, except during Heinrich events.
It is unlikely these anomalously high %Nps levels, corresponding to particularly cold SSTs, reflect a southward movement of the NAPF, when %Nps values further north remain significantly lower. The most likely explanation for the anomalies is that there is colder water during Marine Isotope Stage (MIS) 2 west of France, with this extending closer to the coast during the last few thousand years of MIS 2. This hypothesis is consistent with both reconstructions of North Atlantic SST and some Palaeoclimate Modelling Project 2 (PMIP2) simulations, while not conflicting with the modelling result here of long-term retention of a region of weak flow and divergence off western Ireland. The GLAMAP 2000 project (Pflaumann et al., 2003;Sarnthein et al., 2003) reconstructed both summer and winter SSTs over the whole Atlantic using a large number of sites, and transfer functions, in part, using %Nps. They found a distinct summer cold core area in the central North Atlantic, off southern Ireland and western France, with locally more sea ice in the winter (Figs. 7 and 8 in Pflaumann et al. (2003)), with suspected flows separating around this region, as is consistent with the model circulation here (Fig. 4). Similarly, some, but not all, of the PMIP2 simulations of Li et al. (2010) show cold water, and sea ice spreading eastward from the Newfoundland/Labrador Sea area, leading to similar cooling south of Britain, in water not directly connected to the seas off northwestern Ireland and Scotland. The model used here does not show this cold water tongue, as is also the case for half of the PMIP2 simulations, but does support Pflaumann et al.'s (2003) hypothesis of a weak, diverging flow west of Ireland. In support of this are also the findings of both the current modelling and Scourse et al. (2009), that there is a divide in sourcing IRD throughout non-Heinrich periods of MIS 2 off the Porcupine Bank.
The answer to Question 5 is therefore complex. Changes in climate over time will impact the ocean properties affecting iceberg travel and longevity. However, our current knowledge of the detail of multi-millennial-scale change in these properties, even in the North Atlantic, from both observational and modelling perspectives, is sufficiently limited for it to be problematic to make links between changes in climate, ice sheet and ocean confidently. Despite this, the results described here generate new hypotheses for field and model testing.

Conclusions and consequences
The overall conclusions of the IRD perspective on the peak and deglaciation phases of the last BIISthat IRD is sourced fairly locally over 31-16 ka BP, except during some Heinrich events in some core locations, and that the likelihood of the mean NAPF separating northward-from southward-moving currents remained around the Porcupine Bankhave been shown to agree with the model simulations. These are less supportive of oscillation of this NAPF boundary during MIS 2 and the end of MIS 3, as they show little movement of this boundary during 'normal' glacial conditions across the study period. Earlier, a possible resolution of the %Nps inconsistencies leading to that Polar Front hypothesis is given, through the inferred presence of a cold tongue extending across the Atlantic south of Britain. Nevertheless, during Heinrich events, which is when Scourse et al., 2009 found %Nps to be greatest at all sites, colder conditions are modelled off western Europe (Levine and Bigg, 2008;Bigg et al., 2011), consistent with both the presence of IRD from the LIS and high levels of Nps. Future work using a fully coupled climate model over this time period would help resolve questions about the stability of the NAPF position.
The details of the IRD-model comparison, and the answers to some of the questions posed earlier, highlight a number of points for future consideration. While for much of the time background IRD levels are stable enough for meaningful comparisons to be made with time-slice models, when IRD peaks occur they tend to be short-lived. Even 1000-year time slices proved difficult for models to accurately reproduce the timing of IRD peaks, with often a 1-2000-year mismatch. Obtaining excellent age control of marine cores is vital in such comparisons. Another question the comparison presents is how ice sheet change is reflected in iceberg, and hence IRD, flux. The case of northwestern Scotland presented this question most clearly. Conceptual models based on ice sheet growth lead to maximum model flux at peak growth, but both Rosemary Bank and Barra Fan suggest peak iceberg flux during the growth phase, and to a greater degree, during the rapid retreat phase of the local continental shelf margin phase of the BIIS. The climate model was able to capture this dual signal at Barra Fan, but not at Rosemary Bank, possibly because the Barra Fan signal was a convolution of signals from the Hebridean and Donegal Ice Streams (Fig. 6). The IRD record suggests iceberg flux was linked to the speed of ice expansion and retreat. This goes some way to clarifying the question on the indicative meaning of IRD in terms of ice advance or retreat raised in Scourse et al. (2009) and repeated here as Question 3.
The question over the extent of the ice sheet over Porcupine Bank (Peters et al., 2016) was not resolved by either analysis. The main long distance iceberg travel from calving in this area is to the south, meaning that neither the model nor the IRD record solves this problem. The model suggests calving to the south is the dominant source region for these far-flung Porcupine Bank icebergs, with most calved from any northern source melting, or grounding, locally. Further resolution of this question will need to rely on sea floor morphological evidence.
Much of the British-Irish margin was well reconstructed by the model, with iceberg origins compatible with the IRD signal. Where the largest discrepancy occurred was to the south, at the Armorican Fan. This highlighted that during rapid retreat of an ice sheet into large enclosed seas, such as the Celtic and Irish Sea, modelling needs to represent such coastline changes more rapidly than is commonly done in time-slice experiments. In that case modelled icebergs would have been released both further north and in enclosed seas where there would be great scope for grounding and melting preventing icebergs from joining the off-shelf southward circulation. In this context, development of an additional IRD depositional model within the iceberg model may help with interpretation of IRD records (e.g. Death et al., 2006), but such a model would still mostly help with interpretation of the background IRD record rather than periods of extreme change, when rapid changes in calving, current rearrangement and local bathymetry exceed most climate models' capabilities. Nevertheless, overall our novel study has shown that highresolution IRD records can be successfully linked to temporally evolving iceberg trajectory modelling.

Data Accessibility
Modelling data are available on request from the authors. IRD data are available as a mix of data in public repositories that issue datasets with DOIs and some on request from the authors.

Supporting information
Additional supporting information may be found in the online version of this article at the publisher's web-site. Table S1. Mean fluxes for key circulation parameters during a cross-section of the simulations. Means are taken over the last 2 years of the respective simulation (see Fig. S2 for typical variation of some of these). All fluxes are in Sv and air temperature (T) in°C. ST=sub-tropical, SP=sub-polar, SH=Southern Hemisphere. All dates are ka BP. Figure S1. Variation of orbital parameters over the past 50,000 years. Also shown, in the bottom panel, is the impact on daily July insolation at the top of the atmosphere at 65°N. The period of 31-16 ka BP studied here is marked by the double-headed arrow. Figure S2. Plots of four major large-scale ocean parameters in the North Atlantic region during the transition from the last 10 years of the spin-up (where the orbital parameters are as for 21 ka BP) through 10 years of forcing with orbital parameters for 31 ka BP. Note that by year 5 of the new state the variables are approaching a new annual cycle equilibrium. This is the largest parameter change relative to the spin-up phase. Figure S3. The mean 30 m ocean current over the last two years of a selection of simulations is shown for the NE Atlantic (16, 20, 26 and 30 ka BP). Arrow length is linearly related to speed, with longest~30 cms -1 . The data is shown over the model x and y grids. Figure S4. Bar charts of the IRD source variation over time within 2°of the 5 IRD marine cores, using modelled iceberg mass (in Tg) as a proxy. This simulation also includes the Porcupine and Irish Sea modified release sites, shown in Fig. 7. Note that the dotted line is the total BIIS iceberg flux (in km 3 yr -1 ; see right-hand scale), as reconstructed in Fig. 3. See Figs. 1 and 3 for release sites.