The role of megatides and relative sea level in controlling the deglaciation of the British–Irish and Fennoscandian ice sheets

Key external forcing factors have been proposed to explain the collapse of ice sheets, including atmospheric and ocean temperatures, subglacial topography, relative sea level and tidal amplitudes. For past ice sheets it has not hitherto been possible to separate relative sea level and tidal amplitudes from the other controls to analyse their influence on deglaciation style and rate. Here we isolate the relative sea level and tidal amplitude controls on key ice stream sectors of the last British–Irish and Fennoscandian ice sheets using published glacial isostatic adjustment models, combined with a new and previously published palaeotidal models for the NE Atlantic since the Last Glacial Maximum (22 ka BP). Relative sea level and tidal amplitude data are combined into a sea surface elevation index for each ice stream sector demonstrating that these controls were potentially important drivers of deglaciation in the western British Irish Ice Sheet ice stream sectors. In contrast, the Norwegian Channel Ice Stream was characterized by falling relative sea level and small tidal amplitudes during most of the deglaciation. As these simulations provide a basis for observational field testing we propose a means of identifying the significance of sea level and tidal amplitudes in ice sheet collapse.


Introduction
Ice-ocean interactions are integral to some of the most important feedbacks modulating the global climate system and sea level. The interactions between ocean tides and ice sheets are increasingly recognized as central to these feedbacks, and recent observations and numerical simulations highlight several important mechanisms by which ocean tides influence ice sheet dynamics. Tides have been observed to induce periodic variability in ice stream velocity downstream of the grounding line (Gudmundsson, 2006) and up-glacier, partly by modulating stick-slip motion (Anandakrishnan et al., 2003;Bindschadler et al., 2003;Gudmundsson, 2006Gudmundsson, , 2011Aðalgeirsd ottir et al., 2008;Winberry et al., 2009;King et al., 2010King et al., , 2011. The presence of large tides near the Rutford Ice Stream in the Weddell Sea sector of Antarctica accounts for a 12% increase in ice stream mean flow than would be the case without tides (Rosier et al., 2015). Both elastic (King et al., 2011) and viscoelastic (Gudmundsson, 2011) models have been able to simulate the influence of tides on ice stream velocities, the latter reproducing non-linear interactions. Ice shelves play a critical role in buttressing ice streams, and observations following the recent collapse of ice shelves along the Antarctic Peninsula, such as Larsen B, demonstrate an increase in catchment ice stream velocities following de-buttressing (Rignot et al., 2004;Bamber et al., 2007;Rott et al., 2011). Tides directly influence ice shelf stability through ice fracturing at the grounding line because of vertical displacement and therefore influence iceberg calving (Conway et al., 1999;Schoof, 2007;Goldberg et al., 2009;Vieli and Nick, 2011). Tides also induce mixing beneath ice shelves that modulates basal melt rates Mueller et al., 2012). Rosier et al. (2014) present tidal amplitude simulations for major constituents in response to projected changes in extent and/or thinning of the major ice shelves in the Weddell and Ross Sea sectors of Antarctica. The changes in water column thickness and coastal topography induce an increase in the M 2 tidal constituent of over 0.5 m beneath much of the Rønne-Filchner Ice Shelf, highlighting the potential vulnerability of ice shelves to this tidally induced positive feedback. Increases in tidal amplitudes therefore promote increased ice stream flow rates, promote iceberg calving and destabilize ice shelves, and the de-buttressing that results from ice shelf collapse promotes ice stream acceleration. Tides are therefore integral to a series of positive feedbacks that promote ice sheet collapse with implications for global sea level (Hemming, 2004;Fl€ uckiger et al., 2006).
Tide-ice shelf interactions have been invoked to explain the catastrophic ice sheet collapses integral to the Heinrich events that characterized the glacial North Atlantic (Arbic et al., 2004a(Arbic et al., ,b, 2008 and for rapid deglaciation of the Weddell Sea sector following the Last Glacial Maximum (LGM) (Wilmes and Green, 2014). Several independent global palaeotidal model simulations are now available (Egbert et al., 2004;Uehara et al., 2006;Arbic et al., 2008;Griffiths and Peltier, 2009;Green, 2010) that indicate that the glacial North Atlantic, with ice-equivalent sea levels 120-130 m lower than present (Peltier and Fairbanks, 2006), was characterized by megatidal amplitudes (amplitude > 5 m, tidal range > 10 m) and increased tidal dissipation rates (Griffiths and Peltier, 2009;Wilmes and Green, 2014). These simulations replicate megatides in the same locations around the margins of the mid-and high-latitude glacial North Atlantic, forming a 'horseshoe'-like pattern extending from the Labrador coast-Hudson Strait, across Davis Strait, southern Greenland, Denmark Straits, Iceland, Faeroe-Shetland Channel, and southwards along the European Atlantic margin as far as North Africa (Fig. 1). The glacial North Atlantic was subject to megatides because the ocean basin À devoid of fringing shelf seas À was close to resonance with the M 2 tidal constituent (Egbert et al., 2004;Uehara et al., 2006;Arbic et al., 2008;Griffiths and Peltier, 2009;Green, 2010;Wilmes and Green, 2014). The deglaciation of Hudson Strait at around 8 ka BP was a key event in moving the basin away from resonance with the semi-diurnal tide; North Atlantic tides after 8 ka BP are significantly damped and similar to present (Uehara et al., 2006;Wilmes and Green, 2014). These results support the damped harmonic oscillator hypothesis for North Atlantic tides (Egbert et al., 2004;Green, 2010). Griffiths and Peltier (2009) and Wilmes and Green (2014) identify the LGM grounding line position in Antarctica as an important accessory control, with enhancement of the North Atlantic M 2 tide and dissipation rates with a simulated grounded, rather than floating, Antarctic Ice Sheet margin.
The predicted locations of megatides in the glacial North Atlantic adjacent to some of the most important ice stream discharge locations of the LGM ice sheets, including Hudson Strait (Laurentide Ice Sheet), south Greenland, Iceland and the British-Irish Ice Sheet. Notably, however, the Atlantic margins of the Fennoscandian Ice Sheet, lying to the north and east of the Faroe-Shetland Channel, were micro-or mesotidal during the LGM (Fig. 1). Given the significance of tides in controlling grounding line and fringing ice shelf stability, iceberg calving rates and ice stream dynamics, it is reasonable to infer that the model-replicated patterns of LGM megatides and tidal amplitude change during deglaciation modulated rates of ice sheet disintegration and collapse in different locations and through time. In a recent case study, Bayesian modelling of geochronological data constraining the retreat of the Irish Sea Ice Stream revealed changes in retreat rate during different phases of deglaciation (Chiverrell et al., 2013); rates of 500 km a À1 were attained soon after maximum ice extent at 24 ka BP (Scourse et al., 2009a), but slowing and oscillating thereafter as the ice stream retreated northwards into the Irish Sea Basin. Ice sheet and ice stream retreat rates, such as those of the Irish Sea Ice Stream, are a function of a combination of external variables, including atmospheric heat flux and moisture supply, marine climate, topographic constraints (including adverse slopes in marine-based settings), changes in relative sea level and tidal amplitude alongside internal ice dynamics. Until recently it has not been possible to isolate the tidal controls on ice sheet retreat rates, but high-spatial-resolution palaeotidal model simulations are now available that identify the time-dependent changes in tidal amplitudes for known ice stream discharge locations. These simulations generate hypotheses for future field testing by highlighting possible tidally induced accelerations/decelerations in retreat rate along ice stream axes that can be tested by the acquisition of geochronological data ( 14 C, terrestrial cosmogenic nuclide rock-exposure dating, optically stimulated luminescence) from key localities.
The aim of this paper is to present the results of two highresolution palaeotidal model simulations for the NW European shelf seas from the LGM through the last deglacial cycle alongside published glacial-isostatic adjustment (GIA) predictions of relative sea-level change for key ice stream drainage axes of the British-Irish and Fennoscandian ice sheets. Contrasts in relative sea level and palaeotidal amplitudes between different ice streams through time are identified and hypotheses for future field testing are erected. The presence of megatidal conditions at ice-marginal locations during the last deglaciation is highlighted. There are no analogues for such settings at present on Earth, with the highest maximum tidal range of 8.86 m experienced at the grounding line of the Rutford Ice Stream in the Weddell Sea based on dataassimilating Antarctic tidal modelling (Padman et al., 2008).

Methods
Data from two independent palaeotidal models for the tidal evolution of the NW European shelf seas are presented, denoted here KUTM-Lambeck and ROMS-Bradley (see below for definitions). Both models generate simulations for thousand-year time steps from the LGM to the present (last 22 000 years) and the timescales used in both are based on calibrated radiocarbon ages and so approximate the calendrical timescale. The first simulation is described by Uehara et al. (2006) and in subsequent applications (Scourse et al., 2009b;Mortimer et al., 2014). M 2 palaeotidal amplitudes were  (2018) generated by parameterizing a two-dimensional finite-difference version of the Princeton Ocean Model (Blumberg and Mellor, 1987), also known as the Kyushu University Tidal Model (KUTM; Uehara et al., 2006), with dynamic palaeotopographies. Model resolution is 1/12˚; full specifications are given in Uehara et al. (2006). Vertical changes in overall crustal elevation are accommodated as a function of glacioisostatic deflection by integrating changes in palaeotopography generated by glacio-isostatic (un)loading using the GIA model simulations based on Lambeck (1995;revised in Uehara et al., 2006) as described by Scourse et al. (2009b). The palaeotopographies for each time step combine the modern shelf bathymetry with changes in relative sea level generated by the GIA model output. The bathymetry of the present-day NW European shelf area was prepared by combining datasets from several sources including the UK Proudman Oceanographic Laboratory, British Geological Survey and the US National Geophysical Data Center.
The ROMS-Bradley model (Ward et al., 2016) was developed using the three-dimensional Regional Ocean Modeling System (ROMS) (Shchepetkin and McWilliams, 2005) with a spatial resolution of approximately 1/24˚longitude and 1/40l atitude; full specifications and calibrations of model output with modern observations are provided by Ward et al. (2016). To facilitate comparison with the KUTM-Lambeck simulations only the M 2 results are presented here, although the full palaeotidal simulations also include S 2 and N 2 (Ward et al., 2016). The M 2 principal lunar semi-diurnal constituent is dominant across the shelf and describes the first-order variability in tidal amplitudes through time; it therefore captures the time-dependent variation in tidal ranges central to the aims of this study. The Lambeck (1995;Uehara et al., 2006) GIA model implements an ice model that is now known on the basis of recent glacial geological and geomorphological evidence to significantly underestimate the magnitude of the last British-Irish Ice Sheet. The more recent Bradley GIA model (Bradley et al., 2011) was therefore used in the ROMS-Bradley simulations to provide dynamic palaeotopography. This model integrates an ice model consistent with the then available evidence on the timing, extent and vertical thickness of the British-Irish Ice Sheet. The pre-LGM glacial-deglacial history is not included within the Bradley GIA model and hence simulations are only considered for the period 21 ka BP to present as dynamic palaeotopography is not available for earlier time steps.
Both sets of simulations use output from the same global ocean tidal model, based on the Princeton Ocean Model (Uehara et al., 2006), to force the shelf tides. The modern bathymetry for the global model is based on ETOPO2 (Smith and Sandwell, 1997) for the area north of 72˚S and depth deeper than 1000 m and GEBCO Digital Atlas (Centennial Edition) for the remaining area. The palaeobathymetries used in the global model were prepared as described by Uehara et al. (2006) and integrate the global GIA simulations based on the ICE-5G (VM2) model version 1.1 downloaded in December 2006 (Peltier, 2004). We do not use ICE-5G (VM2) for the NW European shelf GIA because the Lambeck (1995) and Bradley et al. (2011) GIA simulations are at significantly higher spatial resolution. Ward et al. (2016) provide a detailed comparison of the KUTM-Lambeck and ROMS-Bradley models and conclude that ROMS-Bradley provides a better fit to all the available observational constraints of the input hydrodynamic and GIA models. Further comparisons are presented below, but the analysis focuses on the ROMS-Bradley output.
Time-dependent simulations for relative sea level and palaeotidal amplitude are presented from both models for 'sampling points' along the axes of the major ice streams draining the last British-Irish Ice Sheet and for the Norwegian Channel Ice Stream of the Fennoscandian Ice Sheet. The locations and deglacial history of these ice streams are based on glacial geological observations (e.g. megascale glacial lineations, trough mouth fan locations) and ice sheet model evidence compiled by the BRITICE project team and presented by Clark et al. (2012) in their fig. 17. The reconstructed timings of ice stream margins for the respective sampling points during deglaciation are also based on Clark et al. (2012).

Results
The spatial patterns of relative sea level, tidal amplitude and ice sheet changes for both KUTM-Lambeck and ROMS-Bradley models are presented in Fig. 2. The relative sea level and palaeotidal amplitude sampling points are presented in Fig. 3 superimposed on the British-Irish Ice Sheet deglaciation Scenario 1 presented by Clark et al. (2012). This enables the timing of ice-marginal positions for the ice streams, based on the evidence available to Clark et al. (2012), to be compared with the simulation outputs for both relative sea level and palaeotidal amplitude for the respective locations of the sampling points. Ward et al. (2016) demonstrate that the key input variable influencing the palaeotidal evolution for any location À assuming open ocean tidal forcing is already applied (Uehara et al., 2006) À is the GIA boundary condition. Ward et al. (2016) compare ROMS and KUTM hydrodynamic models forced using both the Lambeck and the Bradley GIA models, and it is clear that the palaeotidal response is similar for any individual GIA model using different hydrodynamic models. Given that the ice model used in the Bradley GIA simulations most closely approximates the deglacial dynamics of the British-Irish Ice Sheet, the outputs from ROMS-Bradley form the focus of the analysis in the discussion. The combined relative sea-level/palaeotidal amplitude plots (Fig. 4) show relative sea-level change output from the Bradley GIA model with tidal amplitudes output from ROMS against elevation and amplitude axes through time. The plots show the elevation of present mean sea level (0 m) and the depth of the current sea bed for the selected locations, enabling emergence to be identified; ice cover derived from ice models included in the Bradley GIA input is also indicated.
To facilitate comparisons of sites within and between ice streams, and to identify/predict when the combined forcing of relative sea level and tidal amplitude might have had a significant positive or negative impact on rates of deglaciation, a 'sea surface elevation index' (SSEI) has been calculated. This index is based on the assumption or inference À based on the observations of multiple interacting processes detailed above À that both rapid relative sea-level rise and high tidal amplitudes will promote rapid marine-based deglaciation, and conversely, that relative sea-level fall and microtidal range will mitigate rapid deglaciation. The SSEI is weighted such that the direction (rise, fall) and rate of relative sea-level change is accommodated, but that the high frequency tidal amplitude is probably more important in controlling calving/melting than the relatively long-term relative sea-level change. The index is calculated as: High tidal amplitudes are therefore able to mitigate against a low rate of relative sea-level fall. Note that this unitless index only considers relative sea-level and tidal amplitude change as forcings of deglaciation. The aim is to generate hypotheses to identify/predict where and when relative sea level and tidal amplitudes combine to promote deglaciation; if these predictions are falsified by the field evidence then this implies that other forcings of deglaciation, i.e ocean and atmospheric heating, topographic and internal ice dynamic controls, are probably dominant for the space-time location  (Uehara et al., 2006) and (B) ROMS-Bradley (Ward et al., 2016). The colour scale denotes tidal amplitudes (half full tidal range) from zero (amphidromic) in blue to 4 m in red. Time steps are shown in white text on each panel (a BP, 0 is the present day). Note the increase in megatidal amplitudes along the European margin at the LGM. Ice sheet extent [based on ice models used by Lambeck (1995) and Bradley et al. (2011)] are shown in white and land in grey, with the present-day coastline given for reference (black line). In general, the western ice streams draining the British-Irish Ice Sheet were characterized by much higher tidal amplitudes than those draining eastwards into the North Sea or from the Fennoscandian Ice Sheet. The contrast in the tidal evolution between the northern and southern sectors of the Celtic and Irish seas -already indicated in the KUTM-Lambeck output -is reinforced by the ROMS-Bradley simulations. Scourse et al. (2009b) postulate that the megascale linear sediment ridges of the Celtic Sea represent tidally remobilized sediment sourced from the receding Irish Sea Ice Stream marine margin and the Fleuve Manche because of high tidally induced bed stresses during deglaciation. The contrast in deglacial tidal dynamics across the sediment ridge field described here adds a complication to this hypothesis, implying that either these ridges are not primarily of tidal origin, or that only those in the southern sector are of dominantly tidal origin. If the latter is the case then this suggests that the northern sector of the ridge field formed through non-tidal processes but with an equifinal resulting morphology. This seems unlikely. Alternatively, both of the palaeotidal models might underestimate tidally induced bed stresses in the northern Celtic Sea sector; this could result from GIA uncertainty in this region possibly as a function of ice model inaccuracies. More observational data on the origin and timing of the formation of these sediment ridges are required relating, in particular, to any morphostratigraphic contrasts between the northern and southern sectors of the ridge field.
The most notable feature of the SSEI outputs (Fig. 5) for each sector are the high values generated for the periods between 16-12 ka BP and 11-9 ka BP for all sectors. The latter period will only be relevant for any late Younger Dryas ice remaining with a marine terminus in the innermost parts of sea lochs in Scotland, but the earlier phase, which immediately follows Heinrich Event 1, may well have contributed to rapid deglaciation in the northern Irish Sea north of the Heinrich Event 1 limit (cf. Clark et al., 2012;Chiverrell et al., 2013) in the Barra-Donegal sector and in the inner parts of the Minch sector. If analysis of geochronological data ultimately demonstrates a coeval acceleration of  (Ward et al., 2016). Asterisks denote ice front recessional positions as presented by Clark et al. (2012). Note that the scale for the tidal amplitudes (red lines, tidal amplitude above and below mean sea level) is given in the top panel. Note the difference between the Norwegian Channel Ice Stream low deglacial tidal amplitudes compared with the high tidal amplitudes simulated for the western British-Irish Ice Sheet ice streams. deglaciation between 15.5 and 12.5 ka BP then these simulations suggest this might be attributable to the combined effects of relative sea level and palaeotidal amplitude changes.
The SSEI data enable predictions of phases of rapid deglaciation based on the combined effects of rapidly rising relative sea level combined with high tidal amplitudes. The highest SSEI value of the entire dataset is 14.2 for the Galway sector at 14.5 ka BP at location 3. This is close to the known position of the ice front based on existing data and, along with high values for outer to mid-shelf locations in the Barra-Donegal and Celtic Sea sectors, predicts  (2018) forcing of rapid deglaciation at this time. The Minch sector shows somewhat antiphase behaviour, with a high index at 13.5 ka BP as the ice front reaches the inner parts of the Minch and lower values earlier in deglaciation. This represents another hypothesis to test against new observational data. The very high M 2 amplitudes simulated for the Minch sector are to some extent mitigated by relative sealevel fall driven by glacio-isostatic uplift. The Norwegian Channel Ice Stream has some of the lowest SSEI values of the dataset and in any case the highest values in this sector do not coincide with the known position of the ice front; here, rapid deglacial relative sea-level fall along with low simulated tidal amplitudes combine, suggesting low influence on deglaciation from these forcings.
The reduction in open North Atlantic amplitudes following the deglaciation and consequent flooding of Hudson Strait around 8 ka BP suggests that an analogous reduction in amplitudes might be expected during Heinrich events, since the enhanced iceberg flux that defines these often originate from this sector of the Laurentide Ice Sheet (Hemming, 2004). Such an amplitude dampening is not observed in the M 2 simulations presented here because ICE-5G (VM2) is consistent with the available geological evidence in incorporating no major ice-marginal recession across Hudson Strait during Heinrich Event 1. Establishing ice-marginal locations and the extent of ungrounding across Hudson Strait during Heinrich events is therefore an important geological priority because small changes in this sensitive region probably have implications for the palaeotidal evolution of the wider North Atlantic.

Conclusions
Multiple palaeotidal amplitude simulations demonstrate that the glacial North Atlantic was megatidal. High-spatialresolution palaeotidal simulations for NW Europe enable, for the first time, the influence of changing tidal amplitudes within specific ice stream sectors draining the British-Irish and Fennoscandian ice sheets to be analysed through time.
The simulations predict high tidal amplitudes during deglaciation for the British-Irish Ice Sheet ice streams draining westwards into the North Atlantic, but the results are spatially heterogeneous. The Celtic Sea, Galway Bay and Barra-Donegal sectors are characterized by megatidal amplitudes during either early or mid-deglaciation, whereas the Minch Ice Stream is characterized by megatidal amplitudes during the final stages of deglaciation. Since the palaeotidal simulations use relative sea level as an input variable based on dynamic topography from GIA models, relative sea level can also be isolated along with tidal amplitude as a potential forcing of deglaciation. Relative sea level and tidal amplitude have been combined into an index of 'sea surface elevation' (SSEI) that enables the influence of these combined forcings to be analysed for each ice stream sector through time. This index highlights a marked difference between the western British-Irish Ice Stream ice stream sectors, in which combined relative sea level and tidal amplitudes were potentially important drivers of deglaciation, and the Norwegian Channel Ice Stream where relative sea level was falling and tidal amplitudes were small during most of the deglaciation. These simulations provide hypotheses for future field testing via the acquisition of high-quality geochronological data constraining deglaciation.