Reconstruction of the palaeo‐sea level of Britain and Ireland arising from empirical constraints of ice extent: implications for regional sea level forecasts and North American ice sheet volume

Reconstructions of palaeo‐sea level are vital for predicting future sea level change and constraining palaeo‐ice sheet reconstructions, as well as being useful for a wide array of applications across Quaternary Science. Previous reconstructions of the palaeo‐sea level of Britain and Ireland relied on a circular tuning of glacio‐isostatic models: input ice sheet thicknesses and extents were iteratively altered to fit relative sea level data. Here we break that circularity by utilizing new data from the BRITICE‐CHRONO project, which constrains the position of the British–Irish ice sheet margin through time, and we compare derived glacio‐isostatic modelling to the rich relative sea level record. We test a combination of plausible ice thickness scenarios which account for the uncertainty of ice margin position over the North Sea, demonstrating the region where regional sea level data could distinguish between different glaciation scenarios. Our optimal reconstruction is then combined with several global‐scale reconstructions. As the signal of the British–Irish Ice Sheet is constrained, we demonstrate how the relative sea level record of Britain and Ireland can be used to test reconstructions of far‐field ice sheets (e.g. Antarctica, Eurasia and the Laurentide). The derived palaeo‐topography data are likely to be useful for multiple disciplines. Finally, our improved method of sea level reconstruction impacts predictions of contemporary vertical land motion.


Introduction
The waxing and waning of ice sheets raised and lowered sea level throughout the Quaternary (Chappell et al., 1996;Lea et al., 2002;Spratt & Lisiecki, 2016). These past changes have an ongoing impact, as the Earth continues to deform in response to deglaciation after the Last Glacial Maximum (~26-21 ka BP) (Schumacher et al., 2018). The fingerprint of sea level change induced by ice sheet growth and decay has a local and a far-field signal (Clark, 1976;Milne et al., 2009), recorded in fossils and sediments which act as sea level index points (SLIPs) (Shennan, 1986;Khan et al., 2019). Reconstructing palaeo-ice sheets and identifying their sea level fingerprint enables us to study the long-term role of ice sheets in the Earth system (Stokes et al., 2015), and is crucial for understanding present-day sea level changes and vertical land motion (VLM) (Whitehouse et al., 2012;Mazzotti et al., 2011;Mey et al., 2016).
Broadly, three approaches to reconstructing the extent and dynamics of palaeo-ice sheets exist: (i) empirical reconstructions, whereby the margin position through time is inferred from landforms, sediments and the timing of ice-free conditions from dated materials (e.g. Bentley et al., 2014;Clark et al., 2012;Hughes et al., 2016;Dalton et al., 2020); (ii) ice sheet modelling, whereby the dynamics of palaeo-ice sheets are simulated using a physics-based numerical model driven by climatic fluctuations (e.g. Tarasov & Peltier, 2004;Hubbard et al., 2009;Patton et al., 2017); and (iii) glacio-isostatic adjustment (GIA) modelling, an iterative approach whereby different ice-sheet reconstructions produced from (i) or (ii) are input into a GIA model and adapted to maximize fit to relative sea level (RSL) data (e.g. Bradley et al., 2011;Lambeck et al., 1998;Peltier et al., 2015). Each of these approaches has its merits and pitfalls. For instance, though the empirical approach adheres to data on the ice margin position through time, only a two-dimensional line-on-map reconstruction is provided. The empirical field evidence of vertical ice sheet dimensions (e.g. from trimlines) are spatially limited and open to contrasting interpretations (Ballantyne, 2010;Ballantyne, 2015). Conversely, the ice-sheet modelling simulations provide three-dimensional reconstructions but do not necessarily adhere to, and may vary widely from, data on ice margin extent. In GIA modelling, there is a causality problem: icesheet thickness reconstructions are refined to improve the model fit to RSL data, in a manner that might not be consistent with available geomorphological evidence or glaciological principles. Thus, the integration of these approaches has been a long-term goal of the field (Andrews, 1982;Stokes et al., 2015), and efforts to better integrate data across the approaches are increasingly common (Tarasov & Peltier, 2004;Briggs et al., 2014;Patton et al., 2016;Ely et al., 2021;Gandy et al., 2021).
A number of studies have combined methods of ice sheet reconstruction (e.g. Kuchar et al., 2012;Auriac et al., 2016;Gowan et al., 2021). For the British-Irish Ice Sheet (BIIS), the results of Kuchar et al. (2012), who used the model simulations of Hubbard et al. (2009) as an input for a GIA model, supported the reinterpretation of trimlines across Scotland as a minimum constraint on ice surface altitude. However, the approach of Kuchar et al. (2012) did not directly utilize an empirical reconstruction. A fully hybrid approach was presented by Gowan et al. (2021). They used a reduced icephysics model (Gowan et al., 2016) to generate an input ice thickness reconstruction for a GIA model. The simplified model reduces the number of input parameters, compared to more complex ice-sheet models (e.g. Hubbard et al., 2009), requiring only two main inputs: a map of basal shear stress which alters the ice-surface profile, and empirically derived ice-margin positions which prescribe the extent of ice through time. Since the ice margin can be independently constrained by empirical reconstructions (e.g. Hughes et al., 2016;Batchelor et al., 2019), this reduces the circularity inherent in the usual GIA approach whereby the ice thickness model is tuned to RSL constraints.
The last BIIS provides a data-rich environment to test this new hybrid approach to GIA modelling, with a rich RSL database against which GIA models can be constrained (Shennan et al., 2018). Previous work has achieved good fits between GIA models and RSL data (e.g. Bradley et al., 2011), but used the circular GIA approach. Subsequently, through the BRITICE-CHRONO project , a high density of geochronological data constraining the retreat of the BIIS has been collected (Bradwell et al., 2021a, b;Evans et al., 2021;Chiverrell et al., 2021;Scourse et al., 2021;Ó Cofaigh et al., 2021;Benetti et al., 2021). This has culminated in an empirical reconstruction of ice margin history from 32 to 15 ka BP at 1-ka time steps (Clark et al., 2022a). These ice margin positions supersede those used in previous GIA models (Bradley et al., 2011). Despite the increase in available data, uncertainty in the ice margin position through time remains, especially in comparatively data-poor regions such as the North Sea (Clark et al., 2022a). Furthermore, the thickness of the ice sheet remains poorly constrained, especially given the reinterpretation of trimlines as possible thermal boundaries rather than ice elevation limits (Ballantyne & Stone, 2015). GIA modelling, which can provide a constraint on ice thickness, may be able to discriminate between different reconstructions of North Sea ice cover, and elucidate plausible ice thicknesses. This is the approach we take in this paper.
The position of the BIIS creates an opportunity to derive information about other ice sheets during the last glacial period (Bradley et al., 2011). Across the British Isles, the RSL information recorded in sea level data (index points and limiting data) is influenced by signals from both local [BIIS and Fennoscandian and Barents Sea ice sheets (EuIS)] and nonlocal [North American ice sheet complex (NAISC) and Antarctic (AIS) and Greenland ice sheet (GrIS)] ice sheets (Fig. 1). The relative importance of local and non-local ice sheets varies between sites across the British Isles. Toward the interior of the ice sheet, the RSL signal is equally sensitive to the loading from the local ice sheet and the non-local signal, as typified by the Forth Valley example in Fig. 1(e). Further south, non-local ice sheets dominate the signal. For example, in Hampshire the NAISC was the main contributor to RSL change here (Fig. 1f). Thus, if we can constrain the local signal of the BIIS, we may be able to evaluate reconstructions of other ice sheets.
Here, we test an initial suite of six different empirically constrained local BIIS ice sheet reconstructions, accounting for uncertainties in the duration and timing of ice extents over the North Sea and plausible ranges of ice sheet thickness. Exploration of regional and temporal differences in the mismatch between the original six input reconstructions and the RSL data led to the development of an improved reconstruction by tuning aspects of the reconstruction and map of basal shear stress. These experiments were important for deciding the optimum ice sheet reconstruction over the North Sea that underpinned the recent BRITICE-CHRONO ice sheet reconstruction (Clark et al., 2022a) and for the production of palaeotopographies used in that work, and which are made available for others in a data publication (Clark et al., 2022b). Using this new constrained ice sheet reconstruction enables us to demonstrate the utility of the UK and Irish sea level database for studying the North American and Antarctic ice sheets. As such, we address the following questions:

Local ice sheet reconstructions (BIIS and EuIS)
The thickness of the BIIS and EuIS through time was modelled using ICESHEET1.0 (Gowan et al., 2016). This plastic ice sheet model requires two inputs: basal shear stress and ice margin position. For the basal-shear stress map, we divided the domain of the BIIS and EuIS into five distinct geomorphological categories: (i) palaeo-ice streams; (ii) marine sediments; (iii) thick sediments; (iv) discontinuous sediment; and (v) exposed bedrock (Supporting Information Fig. S1). Changing the basal-shear stress values increased or decreased the basal slipperiness altering the ice surface slope and ice thickness. As such, we used two basal-shear stress maps, using our estimated minimum (thin ice) and maximum (thick ice) estimates of shear stress for the five assigned categories ( Table 1). The thick ice estimate reaches a maximum thickness of~1450 m along the NW region of Scotland (see Fig. 3), whilst in the thin ice estimate this is reduced to 1250 m. These should not be taken as definite limits for the thickness of the BIIS, but an approximate range that is tested here. Basal shear stress remained constant in time throughout our experiments. For the EuIS, the most credible isochrones of Hughes et al. (2016) were used, and these do not change throughout the experiment. The only region of the EuIS that varies in the paper is across the North Sea, where it coalesces with the BIIS. For the BIIS, three reconstructions were developed arising from results from the BRITICE-CHRONO project as they became available and with consultation with the then available literature (see Table S1). These capture the uncertainty of glaciation duration across the North Sea, and are referred to as Early, Middle and Late. Figure 2 graphically represents the differences between these three ice margin reconstructions, and maps of the key time slices are shown in Figure 3. In the Early reconstruction, ice expands across the northern North Sea, resulting in a coalescence of the BIIS and EuIS at 29 ka BP (Sejrup et al., 2009;Fig. 3a) with ice expanding south reaching a maximum at 27 ka. The maximum extent is short-lived, with ice retreating from 24 ka, splitting along a calving bay, following the reconstruction of Bradwell et al. (2008) (Figs. 2 and 3b). In the Middle reconstruction, ice coalesces first in the Southern North Sea, following Clark et al. (2022a), at 27 ka (Fig. 3c). From this coalescence the ice grows slowly across the North Sea, reaching a maximum extent at 23 ka (Clark et al., 2012;Fig. 2). Splitting of the BIIS and EuIS then occurs after 19 ka, with rapid retreat to the onshore sectors (Sejrup et al., 2016). In the Late reconstruction, the BIIS and EuIS expand slowly out from land-terminating margins across the northern and central North Sea (Fig. 3e), but do not fully coalesce until full glacial conditions at 23 ka. The same rapid retreat at 19 ka as in the Middle reconstruction occurs in the Late model (Sejrup et al., 2016). In all three scenarios, the margin position is the same from 18 ka BP onwards.
Using these three ice sheet margin reconstructions (Fig. 2), and the two shear stress input maps (Fig. S1), six initial BIIS and EuIS were generated (Table 2). This initial suite, developed prior to the final reconstruction of Clark et al. (2022a), was used to inform our final 'hybrid' reconstruction. Figure 4 provides an overview of the modelling procedure used to generate this reconstruction. The hybrid reconstruction used the optimum reconstruction from Clark et al. (2022a) as a margin input, and most closely resembles the 'Middle' reconstruction for the North Sea. Though our aim was to keep the margin input and shear stress independent from the RSL comparison, our results showed a discrepancy across Scotland resulting from the spatial pattern of the minimum/maximum estimated shear stress values (see 'How thick was the British-Irish Ice Sheet?' below; Fig. 5). The 'thin' model resulted in an underprediction across NW Scotland, whereas the 'thick' model resulted in an overprediction across Central Scotland. Using these sites, a 'hybrid' basal shear stress map was produced which combined the maximum and minimum patterns (see Fig. S1).

Non-local ice sheet reconstructions
Four different global ice sheet reconstructions were combined with our BIIS and EuIS reconstructions, and input into our GIA model (Tables 2 and 3; Fig. 4). We chose reconstructions which were freely available (early 2021), globally consistent and represented different interpretations of Antarctic (AIS), North American and Greenland Ice Sheet histories (NAISC and GrIS) from 122 ka BP to the present day. (Table 3). The reconstructions are referred to as: (i) Bradley (Bradley et al., 2011); (ii) ICE5G (Peltier, 2004); (iii) ICE6G-AIS (Peltier et al., 2015); and (iv) GLAC1D-NAISC (Tarasov et al., 2012). In ICE6G-AIS or GLAC1D-NAISC only the AIS or NAISC+GrIS respectively have changed relative to the ICE5G reconstruction. These different reconstructions vary in their maximum barystatic sea level (BSL); a large versus small AIS (comparing ICE5G to ICE6G-AIS); timing of the end of AIS deglaciation (Bradley to ICE5G); and style and timing of deglaciation of NAISC (ICE5G vs GLAC1D-NAISC) (see Table 3; Fig. S7). Our aim was not to determine an optimum reconstruction, but rather to investigate if the near-field sea level data can be used like far-field sea level to investigate the NAISC and AIS.
The first two reconstructions (Bradley and ICE5G) were combined with the complete set of seven possible regional reconstructions; the latter two (ICE6G-AIS and GLAC1D-NAISC) were only combined with the final hybrid reconstruction. This led to a total of 16 different ice sheet reconstructions (Table 2).

Relative sea level data
We used the UK and Ireland SLIP database (1541 observations) from Shennan et al. (2018) (henceforth referred to as Shennan, 2018) to quantitatively evaluate the ice sheet reconstructions on a regional scale. The data were separated into seven regions ( Fig. 5; Fig. S2), based on the relationship between the local and non-local signal. For example, for regions where the observed sea level signal is above present day, the local (BIIS and EuIS) signal is larger than the non-local signal (Fig. 1e, Western and Central Scotland) throughout the glacial-interglacial period. Regions where the observed signal is below present (Fig. 1f, SouthWest and SouthEast), the nonlocal signal is larger, and these will be used to investigate the North American and Antarctic ice sheets. Although Shennan, 2018 has 1541 SLIPs, the majority cover the last 10 ka BP, with only West and Central Scotland having SLIP records prior  to 11 ka BP (Fig. S2). We do not attempt in this study to produce a local-scale site-based fit as in Bradley et al. (2011) where the ice sheet reconstruction is tuned to fit all local-scale observations such as limiting data, stratigraphic or sedimentologic records. Therefore, we accept that there may be small local-scale misfits at each specific sea level site.
Glacio-isostatic adjustment model RSL and VLM predictions were calculated using the GIA model, which has been previously used in the study region and applied in other global studies (Yousefi et al., 2018;Love et al., 2016). The model was run using a 1D Earth model     (spherically symmetric, self-gravitating Maxwell body) and a global reconstruction of the Late Pleistocene ice sheets. The effect of changing the latter is the focus of the experimental design of this study and can be separated into local and nonlocal (as detailed above) ice sheet components. When solving the generalized sea level equation (Kendall et al., 2005;Mitrovica & Milne, 2003) perturbations to the Earth's rotation vector, time-varying shoreline migration and sea level change in regions of retreating marine-based ice were considered (Milne & Mitrovica, 1998;Mitrovica et al., 2001).
To optimize the free parameters in the 1D Earth model, lithosphere thickness and upper and lower mantle viscosities, we adopt the approach of Bradley et al. (2011) (Yousefi et al., 2018;Love et al., 2016), whereby we conduct a suite of forward models (114 in total) and calculate the misfit between the observed and modelled sea level using chi squared (χ 2 ): where RSL i p and RSL i obs are the ith predicted and reconstructed RSL respectively, calculated at a given latitude, longitude and time; σ i is the 2 sigma elevational error on the RSL i obs and N is the total number of data points. The initial parameter suite considered lithosphere thicknesses ranging between 71 and 96 km, and upper and lower mantle viscosities between the ranges of 0.1 to 1 × 10 21 and 1 to 50 × 10 21 Pa s respectively. The model was run at 512 spherical harmonics, resulting in añ 35-km resolution on the surface (twice as high as that of Bradley et al., 2011).
Using the entire database of 1541 SLIPs, we identified a set of optimum Earth models which produced the lowest χ 2 .
These parameters were adopted for the regional comparison (see Results and Discusion, and Summary). For each of the seven local sub-regions we compared the minimum χ 2 produced from each of the input ice sheet reconstructions and the suite of Earth models ('How thick was the British-Irish Ice Sheet?'). Note that each region may have a different Earth model which could be different from the optimum (see Tables S2 and S3).

Results and Discussion
How thick was the British-Irish Ice Sheet?
To evaluate the choice of BIIS reconstruction, the RSL predictions must capture both the elevated late Devensian SLIP data across central and west Scotland (Fig. 6) and the location and timing of the highstand during the middle Holocene (central). Previous studies (Bradley et al., 2011;Shennan et al., 2018) have emphasized the complication in generating sufficient local ice loading signal to capture the elevated sea level and not conversely producing an overprediction in the Holocene highstand. In this section we use the five sub-regions across Scotland and central UK and Ireland (Fig. 5) to evaluate the spatial pattern of advance and retreat (Fig. 3) and thickness of the input ice sheet. Additionally, we will compare the RSL predictions at each of the 81 sites [nine shown in Fig. 6 (Fig. S6)]. This subset contains 651 sea level index points, but less than 50% are older than 7 ka. Repeating the analysis that is described below using only older than 7 ka data did not change the conclusions regarding the choice of ice sheet reconstructions (see Table S4 and Fig. S5).
The central, thickest ice dome in all the reconstructions extends along the north-western margin of Scotland (Fig. 3 Figure 6. Comparison of predicted and observed sea level from reconstructions of RSL at nine sites across the study region. Sea level observations: SLIPs (blue crosses), freshwater limiting (green triangles, sea level below); marine limiting (blue triangles, sea level above). RSL predictions are shown for the thick (solid lines), hybrid (solid-black) and thin (dashed) combined with the optimum Earth model (see Table S2) and each of the three different North Sea histories are early (red), middle (grey) and late (brown). (results using the other global models are shown in Fig. S6). The horizontal black line marks an RSL of 0 m. [Color figure can be viewed at wileyonlinelibrary.com] covering the central (purple) and west (orange) Scotland data clusters (Fig. 3). These two regions contain the greatest concentration of older data (Fig. S2) and will be more sensitive to detecting the choice of input ice sheet thickness. The ice thickness transects across this region (Fig. 3g) show the different shape of the ice sheet between these regions. The three regions along the eastern margin (red, pink and brown) are also sensitive to the choice of input ice thickness, but as they predominately contain data <10 ka (Fig. S2) the relative differences will be smaller.
Comparing the misfit within these five sub-regions (Fig. 5) it is evident that there is a pronounced mismatch between the thick and thin reconstructions. Across west Scotland (orange box), the thick model (early-middle-late) results in a substantially lower misfit than the thin model. This is highlighted in Fig. 6, where the thick model captures the elevated SLIP older than 12 ka but the thin model underpredicts by~15 m. The reverse occurs across central Scotland (purple box), where the thick model results in a larger misfit and an~15-m overprediction (Fig. 6).
This western mismatch is replicated along the eastern margin, with the thick model producing the lowest misfit across north Scotland (red box) and the thin model across central region (brown box). The difference between the two reconstructions is of course smaller during the Holocene (<10 ka), with the thick model overpredicting the Holocene highstand by~1 m (Fig. 6). In the east Scotland region (pink box) the large variations between the thick-thin model (with a χ 2 of 280 compared to 60) are due to changes in the choice of input North Sea reconstruction (early-middle-late), which we will discuss in the section below.
To address this spatial mismatch, we created the hybrid reconstruction adopting the revised margin history from Clark et al. (2022a) for the advance and retreat across the North Sea (very close to the middle reconstruction) but with an altered central dome configuration. As the transect in Fig. 3(g) shows, in the hybrid reconstruction, the long central dome of the thick and thin reconstruction is altered to a two-dome structure. The predictions from this hybrid reconstruction are bounded by the thick-thin model (Fig. 6) and capture the elevated sea level in the western and central Scotland regions. We again note that the fit across this region could be improved further by small regional-scale variations in the input ice sheet thickness. However, the aim of this study was to produce an ice sheet reconstruction which is not heavily dictated by the SLIP data. We do not aim to produce a local-scale site-based reconstruction.
For example, one local-scale unresolved issue is the underprediction between 11 and 8 ka BP at sites along the southern coast of England (Fig. 6). The underprediction across Scotland (Arisaig) can be to some extent reduced by changes to the input Earth model and by changes in the choice of input global ice sheet reconstruction ('What can we infer about the North American and Antarctic ice sheets?').

Can we detect the signal of North Sea glaciation?
An objective of this study was to identity regions across the study area which are sensitive to changes in the glaciation across the North Sea. From the regional sea level analysis above, there is only a difference in the χ 2 results for sea level data in NE Scotland when using the thick model and the different styles of North Sea glaciation (Fig. 5). However, over 50% of this dataset is from SLIPs at Tay valley (shown in Fig. 6), whose record is clustered around 8 ka with a variation of~7 m. Repeating this analysis without Tay valley (Fig. S3), there is minimal difference in the χ 2 results, which leads to the question: Can we use the current onshore data to detect the North Sea ice sheet? Using Aberdeen as an example NE onshore site (Fig. 7d) we compared the difference in predicted RSL due to changes in ice sheet thickness (grey-dotted, early-thick cf. early-thin), timing of the onset of glaciation (solid purple, early-thick cf. late-thick) and style of glaciation (red-dotted, middle-thick cf. late-thick). The first thing to note is that the three different styles of glaciation produce pronounced differences at older time steps. However, by 6 ka BP the differences are less than 2 m and there are fewer than five SLIPs. In addition, changing the thickness or style of glaciation produced the same response (both positive difference). Using both the SLIP and limiting data, none of the reconstructions can match both these observations at NE Scotland and Aberdeen (Fig. 7b,c). For example, at Aberdeen (Fig. 7c) only the early-thick, or middle/ late-thin capture the SLIP data, but these reconstructions are ruled out by the limiting data.
As changes to the timing of the onset of glaciation produced by far the largest response in the onshore region (early-thick cf. late-thick) we investigate the magnitude of this signal across the offshore region ( Fig. 7d-f). At 24 ka, in the early-thick reconstruction there is an extensive thick ice sheet across the North Sea compared to a very minor expanse in the late-thick reconstruction. This early loading results in the large +120-m bullseye (Fig. 7d). This signal is reversed by 18 ka, with a double bullseye pattern formed, with the early-thick resulting in an~42-m lower RSL due to the earlier onset of ice retreat. By 8 ka, the contrasting bullseye pattern remains but the magnitude reduces due to the ongoing response of the solid Earth following ice sheet deglaciation. Therefore, as it is not possible to use the current onshore sea level data to constrain the North Sea deglaciation, the answer perhaps lies in the offshore. Maps as described could act as a useful guide for sites to target with the largest sensitivity to the changes in the ice sheet.
What can we infer about the North American and Antarctic ice sheets?
As outlined in the Introduction, RSL across the British Isles is driven by both the local and non-local ice sheets, with the latter dominating the signal at sites across the southern UK and Ireland. In this section we evaluate the choice of input global ice sheet reconstruction, to ascertain if near-field sea level can be used like far-field data to investigate far-field ice sheets (i.e. NAISC and AIS). Previous studies have investigated changes in the deglaciation pattern of the AIS to resolve the misfit between RSL predictions and the UK sea level database Shennan et al., 2002). However, no metric was provided to assess these changes.
The sea level data in the two southern data regions (southwest and east) were grouped into 1-ka bins and the mean residual was calculated (Fig. 8). Note there were too few data younger than 4 ka to group into 1-ka bins. A positive residual implies that the RSL was overpredicted (i.e. higher than the observed), so there was too much melt from AIS or NAISC. A negative residual indicates there was not enough meltwater as the RSL predictions are below the observed. By comparing the mean residual from the four reconstruction we assessed the impact of (i) a large vs. small AIS (ICE5G vs. ICE6G-AIS); (ii) a later end for AIS deglaciation [Bradley (2 ka BP) vs. ICE5G (4 ka BP)]; and (iii) style of NAISC deglaciation (GLAC1D-NAISC vs. ICE5G) (see Fig. S7).
The data in each region were found to be relativity insensitive to changes in the size of the AIS reconstruction, which perhaps is not surprising given how similar the signal Figure 7. (a) The difference in the predicted RSL at Aberdeen between the early-thick and early-thin (grey dotted), middle-thick and late-thick (reddashed), and early-thick and late-thick (purple-solid). Panels (b) and (c) are a comparison of the RSL predictions to observed sea level in NE Scotland and Aberdeen for the seven different regional ice sheet reconstructions; the thick (solid) and thin (dashed) reconstructions combined with the early (red), middle (grey) and late (brown) use a 1D Earth model with a lithosphere thickness of 71 km and an upper and lower mantle viscosity of 5 × 10 20 and 5 × 10 22 Pa s respectively. Blue crosses are SLIPs, green triangles are freshwater-limiting (sea level below); grey triangles are marine-limiting (sea level above). Differences in the predicted RSL between the early-thick and late-thick ice sheet reconstructions across the North Sea are shown in panels (d) for 24 ka, from the two models is over the last 8 ka (see Fig. S7). The differences are larger in the two older time periods (data older 9 ka BP) but still within the size of the data uncertainty. For the southwest region, changes in the global ice sheet reconstructions have very little impact on data 7 ka or younger, with mean residual~± 1.5 m. The Bradley reconstruction has the lowest residual misfit to data 6 ka and younger in the southeast region which may be due to the continuing AIS deglaciation until 2 ka. Although not conclusive this suggests that the data require a later timing for the deglaciation of the AIS, as suggested in earlier studies , Bradley et al., 2011, Brooks et al., 2008. In contrast, it has a notably larger residual in both region's data older than 8 ka, which is probably related to the out-dated NAISC reconstruction. In most of the windows older than 6 ka BP, the GLAC1D-NAISC produced a smallertypically negativeresidual, indicating too little meltwater compared to a positive residual (too much meltwater) from ICE5G. This suggests that with this global ice sheet reconstruction the earlier described misfit between 11 and 8 ka would be reduced. The lower meltwater from the GLAC1D-NAISC reconstruction, which causes the underprediction detected in the southern data, is due partly to a smaller total volume at the glacial maximum (Table 2) and the stillstand in the early Holocene (Fig. S7). These results show the British-Irish sea level data are useful for detecting the other large continental ice sheets (NAISC and AIS). One possible avenue we did not investigate was the contribution of the EuIS (probably the Barents Sea ice sheet) to MWP1a (14.5 ka BP) (Lin et al., 2021). By design our local (BIIS and EuIS) reconstruction was constrained to the EuIS chronology from the DATED-1 project (Hughes et al., 2016) in which the Barents Sea ice sheet's retreat is too late to contribute to MWP1a. To incorporate this possibility, and any other scenarios, requires additional margin constraint data and adaptations of empirical reconstructions. We note that such data are currently sparse in the Barents Sea region.

Implications for present data rates of sea level change and vertical land motion
Quantifying the ongoing background rates of VLM and RSL arising from the legacy of the last glaciation is important for planning land use protections in relation to recent climateinduced sea level rise across the British Isles (Palmer et al., 2018). Using the hybrid reconstruction ('How thick was the British-Irish Ice Sheet?') combined with two different global ice sheet reconstructions  ('What can we infer about the North American and Antarctic ice sheets?') we compared predictions of rates of RSL and VLM across the study region (Fig. 9). Following on from the latter section, the non-local signal will also be contributing to the ongoing sea level and VLM change. The spatial pattern is very similar regardless of the choice of NAISC and AIS (Fig. 9a,b) as this is controlled by the local ice sheet. The different non-local ice sheet reconstructions reduced the maximum rate of sea level fall from −1.15 mm a −1 with the Bradley reconstruction to −0.59 mm a −1 with the GLAC1D-NAISC reconstruction.
Compared to the results from Bradley et al. (2011), there is little change in the maximum predicted rate of sea level change with the updated hybrid reconstruction, −1.1 compared to −1.15 mm a −1 respectively. However, the pattern and more importantly the region of sea level fall has expanded, shifting the location of the zero contour across southern Ireland (Fig. 9b,d). This has important implications for regions such as west and southern Ireland which are hereby predicted to be undergoing sea level fall. Note that these results do not include any 20 th century sea level signal.

The palaeotopography of Britain and Ireland
Our relative sea level predictions can also be used to examine the palaeogeographical changes that have occurred over Britain and Ireland. Figure 10 shows example palaeotopography maps and shoreline positions, demonstrating how now submerged regions were previously subaerially exposed. In a manner similar to previous reconstructions (Brooks et al., 2012), we anticipate that these data will be of wide use for multiple disciplines of Quaternary Science. Indeed, the data have already been utilized in ice sheet modelling experiments (Clark et al., 2022a). A complete dataset of palaeo-topographical maps is available from Clark et al. (2022b).

Summary, conclusions and outlook
Previous reconstructions of palaeo-sea level of Britain and Ireland produced a reasonable fit to relative sea level data (Bradley et al., 2011). However, this fit was based upon an ice sheet reconstruction tuned to fit the sea level data, rather than empirically derived from the ice sheet directly. Here, we have broken this circularity, by using independent evidence for palaeo-ice sheet extent and thickness (Clark et al., 2022a). We find a similar level of fit at the ice-sheet scale to the sea level data with this new approach, but now have greater confidence in our reconstruction: the previous comparable reconstruction (Bradley et al., 2011) fits the data for the wrong reasons; our new work fits with an empirically consistent ice sheet. Using this new set of sea level reconstructions, we find that: • Our new results have implications for improving the method of sea level reconstruction and impact on the predicted contemporary GIA-driven vertical land motion (Palmer et al., 2018;Fig. 9). Compared to the previous reconstructions, the hinge-line (contour of zero vertical motion) moves further south across central and southern Ireland. • Our work provides limits on the potential thickness of the BIIS and divide configuration over the majority of Britain and Ireland. However, we find that the existing onshore sea level data are unable to distinguish between reconstructions with radically different thickness and timing of ice cover over the North Sea. Empirical reconstructions (e.g. Clark et al., 2022a) favour specific North Sea glaciation histories (our hybrid reconstruction) due to independent ice history constraints. The other scenarios are made available here. • Because we have independently constrained and isolated the signal of the BIIS, we utilized the rich dataset of palaeosea level across Britain and Ireland to constrain the timing and magnitude of ice volume changes of the NAISC, finding that the far-field ice sheets from the GLAC1D reconstruction provide the best fit to the older (>6 ka) data across the south coast of the British Isles. Thus, the older data from Britain and Ireland can complement the current RSL records from across North America (Engelhart & Horton, 2012;Engelhart et al., 2015).
Overall, the completion of the BRITICE-CHRONO project (Clark et al., 2022a) means that the BIIS is now the world's most empirically constrained retreating ice sheet. These new results unlock a sea level dataset from across southern Britain and Ireland which contains 873 index points, of which 138 are older than 6 ka and can be utilized to investigate the NAIS. This creates multiple avenues for future work. In this study, we do not conduct site-based tuning, and hence some model-data misfit remains a challenge. The British and Irish sea level data could be used to test different glaciation scenarios of far-field ice sheets (e.g. Lin et al., 2021), whilst generating RSL data to further constrain the glaciation of the North Sea remains a priority. Furthermore, the sea level predictions and related palaeo-topographies generated from this work have uses for other disciplines. Their utility has already been demonstrated in palaeo-ice sheet modelling experiments (Clark et al., 2022a) with potential other uses including migrations of biota (including humans) to the changing palaeo-geography of northwest Europe.
Acknowledgements. The project benefited from the PalGlac team of researchers with funding from the European Research Council (ERC) Advanced Grant to CDC (Grant agreement No. 787263) under the European Union's Horizon 2020 research and innovation programme and which supported S.B. This work was also supported by the Natural Environment Research Council consortium grant BRITICE-CHRONO NE/J009768/1. J.C.E. acknowledges support from an NERC independent fellowship award (NE/R014574/1). We acknowledge Professor Glenn Milne at the University of Ottawa for use of computing resources.

Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Supporting information
Additional supporting information can be found in the online version of this article. Table S1. Development of various time slices of the three North Sea reconstructions. Dark blue marks the timing of the maximum ice sheet extent; grey when North Sea becomes ice free. These discrete time slices are used to run ICESHEET (Gowan et al., 2016). DATED refers to Hughes et al. (2016), BC refers to Clark et al. (2022a). Table S2. Minimum chi-squared misfit for the range of ice sheet reconstruction and associated choice of Earth model parameters. Table S3. Minimum chi-squared misfit for each of the seven regions and the range of ice sheet reconstruction. Table S4. The minimum chi-squared misfit calculated for the subregions shown in Fig. 5 and Fig. S6. Panel (a) is all data and (b) is data older than 7 ka BP. The grey shading highlights for each subregion which model performs the best. Using the older data only increases the misfit for each regional ice sheet reconstruction but does not impact on which model performs the best for each region. Table S5. 1D Earth parameters referred to in the paper with reference name, lithosphere thickness, and upper and lower mantle viscosity. Figure S1. Input basal shear stress values and map used for the Gowan et al. (2016) plastic ice sheet model. Figure S2. Frequency distribution of number of SLIP data points in the seven regions, in 1-ka windows. For each region, a graphical example of the style/characteristic of the sea level history is given, which was used as a basis for separating the data. The dashed dark pink lines for East Scotland are the product of removing the Tay valley dataset. Figure S3. Visualization showing the varying regional misfits of the seven local ice sheet reconstructions to the relative sea level database. For each sub-region the minimum chi-squared misfit is calculated using the seven regional ice sheet reconstructions combined with the Bradley global model. The colour circles are the sea level index locations used in each subregion (box with matching colour); white triangles are sea level sites with only limiting points. The stars are using the earlier published ice sheet reconstructions (Bradley et al., 2011). Plot (a) (LHS) is using the optimum Earth model calculated for the entire sea level database (Table S2); plot (b) (RHS) the results are for the optimum Earth model for each region (Table S3) as shown in Fig. 5, but the Tay valley data are removed for the NE Scotland sub-region. When removing the Tay valley dataset, there is no variation in the chi-squared misfit across NE Scotland (pink box). Figure S4. Visualization showing the varying regional misfits of the seven local ice sheet reconstructions to the relative sea level database. For each sub-region the minimum chi-squared misfit is calculated using the seven regional ice sheet reconstructions combined with the ICE5G global ice sheet reconstruction. Additionally, the results are shown using the GLAC1D-NAISC and ICE6G-AIS global ice sheet reconstructions combined with the hybrid regional reconstruction. The colour circles are the sea level index locations used in each subregion (box with matching colour); white triangles are sea level sites with only limiting points. Plot (a) (LHS) is using the optimum Earth model calculated for the entire sea level database (Table S2); plot (b) (RHS) the results are for the optimum Earth model for each region (Table S3). Figure S5. Visualization showing the varying regional misfits of the seven local ice sheet reconstructions to the relative sea level database. For each subregion the minimum chi-squared misfit is calculated using the seven regional ice sheet reconstructions combined with the Bradley global model and data older than 7 ka BP only. The colour circles are the sea level index locations used in each subregion box (with matching colour); white triangles are sea level sites with only limiting points. The stars are using the earlier published ice sheet reconstructions (Bradley et al., 2011). See Fig. S3. Figure S6. Comparison of the predicted and observed sea level from reconstructions of relative sea level (RSL) at nine sites across the study region. Sea level data: index point (blue cross); fresh water limiting data (green triangles) and marine limiting (blue inverted triangle). RSL predictions are shown using the hybrid regional (BIIS+EuIS) reconstruction combined with the four different global ice sheet reconstructions: Bradley (red) ICE5G (grey); GLAC1D-NAISC (brown) and ICE6G-AIS (blue). Each ice sheet reconstruction is combined with the Earth model which produces the minimum chi-squared misfit (see Table S2). Additionally, predictions are shown using the ICE6G-AIS+hybrid+VM5a model (dotted-blue). The horizontal black line marks an RSL of 0 m. Figure S7. Comparison of the barystatic sea level rise from the set of input global ice sheet reconstructions separated into (a) Antarctic ice sheet only, with results from the three different reconstructions: ICE5G (grey line), ICE6G (blue line) and Bradley (redline), and (b) North American ice sheet complex (NAISC) and Greenland (GrIS) with results from the three different reconstructions: ICE5G (grey line), GLAC1D (brown line) and Bradley (red line). Figure S8. Illustration of the differences in the GIA-driven present-day rates (mm a -1 ) of RSL change (a and b) and VLM (c and d). Results are using the same local BIIS and EuIS reconstruction (hybrid) combined with two different global reconstructions: GLAC1D-NAISC (a and c) and Bradley (b and d). The results are using the same 1D Earth model with a lithosphere thickness of 71 km and an upper and lower mantle viscosity of 5 × 10 20 and 1 × 10 20 Pa s respectively. To highlight the sensitivity to changes in the Earth model, the white contour on all plots is using the same Earth model as in Fig. 9, with strong lower mantle viscosity of 3 × 10 22 Pa s.