Sensitivity of palaeotidal models of the northwest European shelf seas to glacial isostatic adjustment since the Last Glacial Maximum

The spatial and temporal distribution of relative sea-level change over the northwest European shelf seas has varied considerably since the Last Glacial Maximum, due to eustatic sea-level rise and a complex isostatic response to deglaciation of both nearand far-field ice sheets. Because of the complex pattern of relative sea level changes, the region is an ideal focus for modelling the impact of significant sea-level change on shelf sea tidal dynamics. Changes in tidal dynamics influence tidal range, the location of tidal mixing fronts, dissipation of tidal energy, shelf sea biogeochemistry and sediment transport pathways. Significant advancements in glacial isostatic adjustment (GIA) modelling of the region have been made in recent years, and earlier palaeotidal models of the northwest European shelf seas were developed using output from less well-constrained GIA models as input to generate palaeobathymetric grids. We use the most up-to-date and well-constrained GIA model for the region as palaeotopographic input for a new high resolution, three-dimensional tidal model (ROMS) of the northwest European shelf seas. With focus on model output for 1 ka time slices from the Last Glacial Maximum (taken as being 21 ka BP) to present day, we demonstrate that spatial and temporal changes in simulated tidal dynamics are very sensitive to relative sea-level distribution. The new high resolution palaeotidal model is considered a significant improvement on previous depth-averaged palaeotidal models, in particular where the outputs are to be used in sediment transport studies, where consideration of the near-bed stress is critical, and for constraining sea level index points. © 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Since the Last Glacial Maximum (LGM), which occurred globally between ca. 26 and 20 ka BP (Clark et al., 2009), the changing ice mass balance between the oceans and land has had a significant role in driving changes in shelf sea tidal dynamics (e.g. Scott and Greenberg, 1983;Belderson et al., 1986;Austin, 1991;Austin and Scourse, 1997;Uehara et al., , 2006Egbert et al., 2004;Arbic et al., 2008;Neill et al., 2010). The maximum volume of the LGM British-Irish Ice Sheet occurred between 27 and 21 ka BP (Chiverrell and Thomas, 2010;Clark et al., 2012;Everest et al., 2013), during which time the majority of the ice sheet was marine-based, extending into the present-day North Sea and the continental shelves of Britain and Ireland, with different sectors of the ice sheet reaching maximal extents at different times (e.g. Sejrup et al., 2005;O Cofaigh and Evans, 2007;Scourse et al., 2009a;Chiverrell and Thomas, 2010;Clark et al., 2012;Praeg et al., 2015).
It is well known that the climate-induced ice/ocean mass redistribution resulting from periodic advance and retreat of ice sheets causes a spatial and temporal response in sea levels (Fairbanks, 1989;Bard et al., 1996;Lambeck and Chappell, 2001;Peltier, 2002b;Church et al., 2013;Lambeck et al., 2014). Lambeck et al. (2002) estimated that global mean sea level rose by~130 m between 20 and 7 ka BP, and relative sea level (RSL) records suggest that between 7 and 3 ka BP, global mean sea level likely rose a further 2e3 m (Lambeck et al., 2004(Lambeck et al., , 2010. Changes in RSL continue after the increase in eustatic sea level, even once ice melting has ceased, since isostatic responses have a very long lag, giving the Earth system 'memory'. The term 'glacial isostatic adjustment' (GIA) is often given to the isostatic deformation of the solid Earth and the geopotential associated with mass distribution (both solid Earth and ice/ocean). The response of the mantle to surface loading occurs on timescales of thousands-to hundreds of thousands of years, and the earth response-function is both depthdependent and spatially variable (Lambeck and Chappell, 2001). Within shelf seas, changes in sea levels such as those that have occurred since the LGM have significant implications for tidal dynamics, including impacting upon tidal amplitudes (elevations and currents), sediment transport and hence upon sediment type and bedform distribution (Belderson et al., 1986;Van Landeghem et al., 2009b;Furze et al., 2014). Sea level changes also influence seasonal stratification (Austin, 1991;Austin and Scourse, 1997;Scourse et al., 2002;Uehara et al., 2006) and hence shelf sea primary production and carbon dioxide uptake.
Because of limited and expensive data and with increased computing/modelling capacity, palaeotidal model simulations have become an increasingly important tool in Earth systems science over the last 15e20 years, both in deep-time applications and in contexts providing strong sequential age control via radiocarbon dating in the period since the LGM. Coastlines and bathymetry (palaeotopography) change in response to RSL changes forced by primarily eustatic, isostatic and tectonic controls. Numerical tidal models implemented with these palaeotopographies simulate past tidal amplitudes and tide-dependant parameters including bed stress, shelf sea stratification, surface tidal currents and the dissipation of tidal energy. There are multiple implications and applications of tidal models which implement these palaeotopographies and simulate past tidal amplitudes and tide-dependant parameters. Predictions of total dissipation of tidal energy in the Earth system, and the distribution of dissipation between the deep and shallow ocean, have implications for the Earth-moon separation and speed of rotation of the Earth (Munk, 1997;Thomas and Sündermann, 1999;Green and Huber, 2013), and the mechanical energy available to drive the Meridional Overturning Circulation (Green et al., 2009). Accurate simulations of palaeotidal amplitude are important for constraining sea-level index points (SLIPs), and therefore correcting glacial isostatic adjustment models of vertical crustal motion (Neill et al., 2010;Scourse, 2013), and assessment of the changing areal extent of the intertidal zone with evolutionary, ecological and archaeological implications (Balbus, 2014;Mortimer et al., 2013). Reconstructions of the distribution of tidal mixing fronts predict organic flux to the seabed, and the distribution of biota, in shelf sea environments (Uehara et al., 2006;Scourse, 2013), and the role of the transient shelf seas in the carbon cycle (Rippeth et al., 2008). Since palaeotidal models generate bed shear stress predictions, they also provide a means of hindcasting shelf sea sediment and habitat dynamics (van der Molen and de Swart, 2001;Scourse et al., 2009b;Van Landeghem et al., 2009a).
A number of numerical modelling studies have suggested that global tides have changed significantly since the LGM (e.g., Belderson et al., 1986;Thomas and Sündermann, 1999;Egbert et al., 2004;Hall and Davies, 2004;Uehara et al., 2006), notably indicating the presence of large tidal ranges in the glacial Arctic Ocean (Griffiths and Peltier, 2008;Griffiths et al., 2009), Labrador Sea (Arbic et al., 2008), North Atlantic (Arbic et al., 2004;Egbert et al., 2004;Uehara et al., 2006;Green, 2010) and in the western Atlantic (Hill et al., 2011). Global ocean tidal models have demonstrated that the glacial North Atlantic was characterised by megatides (tidal ranges >10 m) in some sectors adjacent to major ice stream margins, with feedbacks to ice sheet dynamics and iceberg calving (Wilmes and Green, 2014;Scourse et al., Submitted for publication). A number of attempts have been made to estimate evolving tidal regimes on various continental shelf regions since the LGM, including the northwest European shelf seas (e.g., Austin, 1991;Scourse and Austin, 1995;Uehara et al., 2006;Neill et al., 2010), the East China Sea (e.g., , and the Bay of Fundy and the Gulf of Maine (e.g., Scott and Greenberg, 1983). Uehara et al. (2006) conducted a detailed numerical modelling study into the tidal evolution over the northwest European shelf seas during the last 21 ka. A two-dimensional (2D) finite-difference palaeotidal model was used to test the differences in incorporating the GIA model of Peltier (1994) and a version of the Lambeck (1995) GIA model into the palaeobathymetries. Further, they investigated the impact of setting the boundary conditions of the domain to be fixed to present-day tides ('fixed') or to vary ('open-ocean'), where the latter accounted for changes in ocean tides caused by changes in eustatic sea level and ice-sheet extent. According to Uehara et al. (2006), the timing of the changes in tidal dynamics was sensitive to the local isostatic effects; however the temporal variability in offshelf tides ('open-ocean') since the LGM had the most profound impact. Although there are existing studies on palaeotidal changes of the northwest European shelf seas since the LGM, these are based upon superseded GIA models such as Peltier (1994) and Lambeck (1995), use 2D hydrodynamic models only, and advancements and availability of computing technologies now make simulations of higher spatial resolution more feasible.
Understanding how the shelf seas have responded to past changes in sea level is useful for predicting the likely impacts of future sea-level rise (SLR). Although projections of future (eustatic) SLR contain considerable uncertainties, ranging from 2 to 10 mm y À1 until 2100 (Church et al., 2013), post-glacial sea-level change provides us with an analogue of how these changes are likely to affect shelf sea dynamics. The shallow shelf seas of northwest Europe are well-suited to such studies due to the evolution of the highly energetic tidal regimes of the region, and due to the response to ice loading and unloading. This work builds upon previous studies by developing a palaeotidal model which incorporates dynamic palaeotopography from the most recent GIA model of the region, developed by Bradley et al. (2011). At the time of its development, this new GIA model included all available RSL data for both Britain and Ireland (Milne et al., 2006;Shennan et al., 2006b;Brooks et al., 2008;Bradley et al., 2009), thus including more recent field data regarding the extent of the ice sheets which were not available at the time of the development of the Lambeck (1995) model. Within the Bradley et al. (2011) model, a revised version of the global model of Bassett et al. (2005) was used for the non-local ice-melt signal, which had been tuned to fit RSL data (from China and Malay-Thailand), and the British-Irish Ice Sheet component of this global model had been removed. This global model also incorporated revised deglaciation of the Antarctic Ice Sheet (see Bradley, 2011, for more details). In further improvements, palaeotidal amplitude outputs from Uehara et al. (2006) were used to apply palaeotidal corrections to the locations of five of the SLIPs used to constrain this new GIA model. Fewer SLIPs were used in constraining the earlier model of Lambeck (1995), none of which were corrected for temporal changes in tidal amplitudes (i.e. none had palaeotidal corrections applied). SLIPs are fundamental to constraining GIA models, and tend to form at or around mean high water spring tide (MHWST), rather than at the palaeo mean sea level (Shennan et al., , 2006a. Together with this reference water level, the 'indicative range' is considered, which is the likely range over which the sediment has been deposited, thus defining the 'indicative meaning' of the SLIP. MHWST varies with tidal amplitude, which is influenced by RSL changes (Gehrels et al., 1995;Shennan et al., 2000a;Uehara et al., 2006;Neill et al., 2010); therefore it is incorrect to assume that tidal range is time invariant. Correcting SLIPs for palaeo-changes in MHWST is thus critical for constraining their indicative meaning. Another component of the Bradley et al. (2011) GIA model was an advanced sea-level change model, which solved the generalised sea-level equation (Milne and Mitrovica, 1998;Mitrovica et al., 2001;Mitrovica and Milne, 2003;Kendall et al., 2005) and which incorporated the most up-to-date developments in GIA sea-level modelling. These advancements included the effect of GIA perturbations on the Earth's rotating vector, improved representation of sea-level change in regions of ablating marine-based ice, as well as time-varying shoreline migration .
The new palaeotidal model presented here (which is referred to as ROMS þ Bradley), incorporates the new GIA model of Bradley et al. (2011) and was developed using the Regional Ocean Modelling System (ROMS), which is an advanced open-source 3D hydrodynamic model (Shchepetkin and McWilliams, 2005). This new palaeotidal model is the most comprehensive, high-resolution (double the spatial resolution of previous models), 3D palaeotidal model for the northwest European shelf seas available to date. We use the model to simulate changes in tidal dynamics across the region from 21 ka BP to present day, in 1 ka time slices. The outputs of this new model are compared with those from two additional palaeotidal models for the region (referred to as ROMS þ Lambeck and KUTM þ Lambeck), which incorporate an older GIA model of the region (Lambeck, 1995). The results are then used to revise the changes in MHWST at Arisaig (Scotland), which has an extensive SLIP record.

The northwest European shelf seas
The northwest European shelf is an interesting region for studying the impacts of changing sea levels since the LGM on shelf sea tidal dynamics, because of the complex RSL signal produced by the unloading of both near-and far-field ice sheets. In the region of the British Isles, RSL change since the LGM has largely been a function of crustal rebound due to ice-unloading of the British-Irish Ice Sheet, which had an ice volume equivalent to a global SLR of 2.5 m (Clark et al., 2012). The unloading of more distant ice sheets, including the Fennoscandian and Laurentide ice sheets, has also had an effect (Lambeck, 1995), although the British Irish-and Fennoscandian Ice Sheets were confluent during the LGM (Bradwell et al., 2008). Because of this complex pattern of RSL, and the extensive sea-level index database available, the British Isles have been a focus of many studies of post-LGM GIA (e.g Lambeck, 1993Lambeck, , 1995Lambeck, , 1996Shennan et al., 2000bShennan et al., , 2006bPeltier, 2002b;Bradley et al., 2011).
On the present-day European shelf, the seas most exposed to the Atlantic waters are the Celtic Sea, the Malin Shelf and the northern North Sea (Fig. 1). There are a number of semi-enclosed epicontinental seas on the shelf, including the southern North Sea, the Irish Sea, the English Channel, the Minch and the Sea of Hebrides. In the North Sea, depths are generally between 100 m and 200 m, although the Norwegian Trench -the majority of which has depths around 300 m -reaches a maximum depth of around 700 m adjacent to the Fennoscandian landmass. The Irish Sea is a semienclosed body of water with a north-south trending channel of depth 250 m, bordered to the south by the Celtic Sea, and linked north to the Malin Shelf by the North Channel. The English Channel, which has a depth of less than 50 m, connects the eastern Celtic Sea to the southern North Sea, and is characterised by large linear sandbanks (Dyer and Huntley, 1999).
The tides on the northwest European shelf seas are predominantly semi-diurnal (Pingree and Griffiths, 1978), dominated by the M 2 (principle lunar) and S 2 (principle solar) tidal constituents. Tidal energy propagates from the North Atlantic onto the shelf into the Celtic Sea. Some of the tidal wave propagates into the English Channel and the North Sea, while some passes into the Bristol Channel and the Irish Sea (Pugh, 1987). Part of the semi-diurnal wave is diffracted south and east by northern Scotland, and enters the North Sea. In the Celtic Sea, the Bristol Channel and English Channel, the diurnal tide behaves as a standing wave, but without any tendency to resonance (Pugh, 1987). The tidal range in the Severn Estuary, which reaches a maximum of~12 m, is the second largest in the world, after the Bay of Fundy. It has long been realised that higher-than-average intensity of energy dissipation occurs in the shallow shelf seas around the UK (Flather, 1976;Simpson and Bowers, 1981), where approximately 5e6% of total global tidal dissipation occurs (Egbert and Ray, 2001;Egbert et al., 2004).

Methods
Here, we present a new high-resolution (3D) palaeotidal model, which we refer to as ROMS þ Bradley (Section 3.1). The results presented within this paper focus on these new simulations, but we also compare these ROMS þ Bradley simulations with two other palaeotidal models (Table 1). One is an existing palaeotidal model, which we refer to as KUTM þ Lambeck, and was developed by Uehara et al. (2006) (Section 3.2). A second new model has been newly-developed (ROMS þ Lambeck), for the purpose of intermodel comparison (Section 3.3).

Palaeotidal model 1: ROMS þ Bradley
The new ROMS þ Bradley palaeotidal model was developed using the 3D Regional Ocean Modelling System (ROMS), an opensource, free-surface, terrain-following, primitive equations model (Shchepetkin and McWilliams, 2005). The model was set up with a spatial resolution of 1/24 longitude and variable latitudinal resolution (1/34 e1/56 ) resulting in~2e3 km grid spacing in latitude and longitude, extending from 15 W to 11 E, and from 45 N to 65 N (Fig. 1). The palaeotidal model is forced at the boundaries using surface elevation (Chapman boundary conditions) and the u (eastward) and v (northward) components of depth-averaged tidal current velocities (Flather boundary conditions). For initial present day model calibration and optimisation, boundary conditions are interpolated from the 1/4 resolution global OSU TOPEX/Poseidon Global Inversion Solution 7.2 (Egbert et al., 1994;Egbert and Erofeeva, 2002). The palaeotidal model forcing is derived from the harmonic constants of the global palaeotidal model output of Uehara et al. (2006) (1/2 global resolution,~25 Â 50 km grid spacing), available at 1 ka time slices from 21 ka BP to the present. Forcing the model with the present-day outputs from Uehara et al. (2006) facilitates direct comparison of the palaeotidal model set-up (ROMS þ Bradley) with the outputs from these optimised presentday runs.
The tidal constituents considered in the derivation of the boundary conditions were M 2 , S 2 , along with the N 2 (larger lunar ecliptic semi-diurnal), these being the three dominant tidal constituents of the region. Incorporating the S 2 tidal constituent alongside M 2 produces spring-neap modulation, where in general peak spring tides induce the highest bed shear stress. Further, the N 2 constituent influences the strength of alternate sets of spring and neap tides, a result of the elliptical orbit of the moon (Pingree and Griffiths, 1981). M 4 is the non-linear shallow water overtide of M 2 , and on the northwest European shelf seas, the phase relationship between M 2 and M 4 dominates tidal asymmetry, which is important for net sediment transport on the shelf (Pingree and Griffiths, 1979). The M 4 tidal constituent is also included in the model output of bed shear stress. The near-bed (or bottom) bed shear stress is automatically set within ROMS to be calculated at the mid-depth of the lowest vertical (sigma) layer. The tidal model was set to have ten layers in the sigma coordinate, using the coordinate system of Shchepetkin and McWilliams (2005). The resolution of the layers was increased towards the bed by adjusting the values of the sigma coordinate bottom/surface control parameters in the model runtime options. The option for quadratic bottom drag scheme was implemented, using a bottom drag coefficient of 0.003, which is in line with other ROMS modelling studies of the region (e.g. Hashemi and Neill, 2014;Lewis et al., 2015;Robins et al., 2015).
The present-day bathymetry was derived from GEBCO (General Bathymetric Chart of the Oceans), at a resolution of 1/120 in both latitude and longitude (~1 km grid spacing in our study region). Alternate wetting and drying of land cells was not included, and a minimum water depth of 10 m was applied throughout the domain. Given that the model resolution does not fully resolve intertidal regions, the lack of alternate wetting and drying is considered acceptable at this horizontal resolution, and for this domain size. To account for GIA, it was necessary to add modelled estimations of GIA to the present-day bathymetry. The palaeotopography dataset implemented in the ROMS þ Bradley model was developed by Bradley et al. (2011), and was available at a resolution of 1/12 for 1 ka time slices from 21 ka BP to present. It should be noted that there is no consideration of dynamic bed level change, whereas in reality some areas have experienced considerable post-glacial sediment accumulation. The melting of the British-Irish and Fennoscandian Ice Sheets was a major source of sediment supply to the northwest European shelf seas and played a significant role in the formation of large bedforms (Boulton et al., 1985;Scourse et al., 2009a) and of large sand banks, such as Dogger Bank in the North Sea (Carr et al., 2006).
For each time slice, the lower resolution GIA model output was bi-linearly interpolated onto the high resolution computational grid. This involved calculating interpolated values of the GIA model outputs at each of the computational grid coordinates such that the present-day bathymetry and the GIA model outputs could be combined to generate a palaeobathymetry for each time slice (Fig. 2). By incorporating the GIA model into the bathymetry grid, some regions loaded by former ice sheets were simulated to have been flooded, since the isostatic loading was greater than the global   (1995) Two-dimensional,~1/12 resolution, 0e21 ka BP in 1 ka time slices Uehara et al. (2006) reduction in sea level. Knowledge of the extents of the British-Irish and Fennoscandian Ice Sheets was thus needed, in order to mask out areas that would have been ice covered. A revised version of the global ice sheet model developed by Bassett et al. (2005) was used here, which is consistent with that used by Bradley et al. (2011). Incorporated into this global ice sheet model (1/2 in latitude, 1/4 in longitude, respectively) were recent constraints on regional ice by Shennan et al. (2006b) and Brooks et al. (2008). Since the work of Bassett et al. (2005), further improvements have been made in constraining the extent of the British-Irish and Fennoscandian Ice Sheets (e.g. Bradwell et al., 2008;Chiverrell and Thomas, 2010;Evans and Thomson, 2010;Clark et al., 2012), and the effort is ongoing.
The model was run for 30 days, from which the last 28 days of model output were analysed. The simulated M 2 and S 2 tidal constituents were separated using harmonic analysis (T_TIDE, Pawlowicz et al., 2002), and were compared with harmonic constants from 19 tide gauges within the UK tide gauge network (National Tidal and Sea Level Facility, 2012) (locations shown in Fig. 1). The root mean square error (RMSE) was 14 cm (8% scatter index, which is the RMSE normalised by the mean of the data) in amplitude and 10 (5%) in phase (M 2 ), and 4 cm (7%) in amplitude and 8 (5%) in phase (S 2 ). To validate the tidal current speeds, published data from 22 offshore current meters (see Jones, 1983;Davies and Jones, 1990;Young et al., 2000, for further details) were compared with the simulated depth-averaged current speed at the grid point nearest each offshore current meter location (Fig. 1), which was also analysed using T_TIDE. The RMSEs of the M 2 tidal currents were 5.8 cm s À1 (18%) in amplitude and 18 (12%) in phase, and were 2 cm s À1 (12%) and 14 (7%) in phase for the S 2 tidal currents. Given the low spatial resolution of the tidal forcing, the present-day model was considered to perform well, giving confidence in the hydrodynamic model set-up and the model outputs.

Palaeotidal model 2: KUTM þ Lambeck
Existing palaeotidal model simulations (KUTM þ Lambeck) were used for comparison with the new palaeotidal model developed here. This model, which is described in more detail in Uehara et al. (2006), is 2D, and was developed using a version of the Princeton Oceanographic Model which is referred to here as the Kyushu University Tidal Model (KUTM). The outputs of Uehara et al. (2006) were available at 1 ka BP time slices for M 2 , S 2 , N 2 , K 1 and O 1 tidal constituents. The palaeotopographies for this palaeotidal model were generated using the revised version of the GIA model of the region developed by Lambeck (1995). The forcing at the model boundaries was consistent with the boundary forcing of both ROMS þ Bradley and ROMS þ Lambeck. This palaeotidal model was forced at the boundaries using the global model outputs of Uehara et al. (2006) (M 2 , S 2 and N 2 , K 1 and O 1 ).

Palaeotidal model 3: ROMS þ Lambeck
To compare with the new 3D ROMS þ Bradley and the existing 2D KUTM þ Lambeck models, a second 3D ROMS palaeotidal model was developed (ROMS þ Lambeck). This second new model incorporated the revised version of the GIA model of Lambeck (1995), consistent with the KUTM þ Lambeck model. The same hydrodynamic set-up within ROMS was used in this model as in the ROMS þ Bradley model. For consistency with the shelf-scale palaeotidal model of Uehara et al. (2006), ROMS þ Lambeck was set up with the same minimum water depth, spatial resolution and extent (6 m, 1/12 longitude and 1/16 e1/28 latitude, covering 15 W to 15 E and 45 N to 65 N). As for the ROMS þ

Results
The outputs of the new, 3D palaeotidal model (ROMS þ Bradley) are presented here in detail, as well as comparison with two other palaeotidal models -ROMS þ Lambeck (new) and KUTM þ Lambeck (existing). To assess the model skill, we compare the model set-ups and the resulting differences in evolving tidal dynamics on the shelf by considering the differences in the palaeotopographies generated using the different GIA models (Section 4.1), as well as changes in M 2 tidal amplitude (this being the dominant tidal constituent on the shelf, Section 4.2). In order to assess how changes in tidal currents may have influenced sediment transport patterns with RSL rise since the LGM, we describe and discuss the evolution of simulated bed shear stress vectors on the shelf (Section 4.3). The new model outputs are then used to present a revised sea-level curve for Arisaig, which is the location of the most extensive SLIP record for the British Isles (Shennan et al., 2006b), with data covering over 16 ka (Section 4.4).

Relative sea-level changes
There are significant differences in land/ice and water boundaries in the model palaeotopographies generated using the Bradley et al. (2011) andLambeck (1995) GIA models (referred to as 'Bradley' and 'Lambeck',Figs. 2 and 3). The increased isostatic depression caused by the thicker and more extensive LGM ice sheets in the ice model within the Bradley GIA model (Shennan et al., 2006b;Brooks et al., 2008;Bassett et al., 2005;Bradley et al., 2011) resulted in higher RSL adjacent to the ice sheet loading in comparison to the Lambeck palaeotopographies. At 21 ka BP, the land/ice mask extended considerably farther north over the western shelf in the Bradley simulations. This can be seen in Fig. 3, where between 20 and 11 ka BP, the areal extent of the submerged shelf seas is considerably greater in the Bradley palaeotopographies. The land/ice extended farther south in the Lambeck simulations until 13 ka BP. After 8 ka BP there was no significant difference between the palaeotopographies generated by the two models.
The northern section of the North Sea flooded considerably earlier according to the Bradley GIA model. Fig. 2 suggests that, in the Bradley model, at 21 ka BP much of the northwest European shelf seas was emerged, and the majority of the shelf that was submerged had water depths <50 m. With focus on the newer Bradley GIA model, following initially rapid deglaciation, a significant area of the shelf was submerged by 20 ka BP, with water depths of up to 150 m. By 15 ka BP, the response of the unloading of the Fennoscandian Ice Sheet was apparent, with areas of the southern North Sea re-emerging as land. This response was also apparent in the northern North Sea, where water depths were shallower at 15 ka BP than during earlier times (e.g. 18 ka BP). Between 9 and 8 ka BP, the English Channel opened, and much of the North Sea was submerged. On the whole, the present-day coastline had emerged by 6 ka BP, with little eustatic SLR occurring during the last several thousand years.

M 2 tidal elevation amplitudes
In all three palaeotidal models, simulated tidal elevation amplitudes and phases over the northwest European shelf seas have changed significantly since 21 ka BP. The focus here is on the outputs from the new ROMS þ Bradley simulations, although key differences in comparison to other palaeotidal model simulations are highlighted. The significantly smaller areal extent of the shelf seas at this time will have also contributed to reduced dissipation of tidal energy on shelf seas in comparison to the present day (Uehara et al., 2006). Changes in the interaction between incident and reflected tidal waves, which is influenced by water depths, topography and frictional effects, can cause migration of amphidromic points (a point of zero tidal amplitude) (Hendershott and Speranza, 1971). Thus, with changing sea levels, the location of amphidromic points can be affected, which in turn influences tidal amplitudes.
The changes in the elevation amplitude and phase of the M 2 tidal constituent output from ROMS þ Bradley, for a number of key Fig. 3. Differences between the land/ice and sea boundaries (i.e. masks) for model palaeotopographies generated using the Lambeck (additional mask in orange fill) and Bradley (additional mask in blue fill) GIA models. Orange line: the land/ice mask used by Uehara et al. (2006) and in the ROMS þ Lambeck simulations, with a minimum water depth of 6 m and interpolated onto the (higher resolution) ROMS computational grid. Blue line: the land/ice masks for the ROMS þ Bradley simulations, with a minimum model water depth of 10 m. The shaded grey regions are the present-day coastline/land. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) time slices, are shown in Fig. 4. Very high M 2 amplitudes were observed on the margin of the eastern North Atlantic during deglaciation, due to the North Atlantic tides being closer to resonance during this time (Egbert et al., 2004), consistent with the findings of Egbert et al. (2004) and Uehara et al. (2006). Due to large North Atlantic tides, there were also areas of high M 2 amplitude along the western Irish coast, which are simulated to have peaked in excess of 3.5 m at 19 ka BP, and remained >2.5 m until 12 ka BP, and between 2 m and 2.5 m until 10 ka BP. Uehara et al. (2006) found that the M 2 amplitude in this area was slightly higher (>2.5 m) until 10 ka BP. During deglaciation, there were also high M 2 amplitudes along the north-western extent of the palaeocoastlines on the shelf, along the west coast of Scotland, again in response to the large tidal amplitudes in the North Atlantic. From 21 to 11 ka BP these were >2.5 m and were particularly high (>3 m) and migrated north between 16 and 12 ka BP, as the shelf in this region flooded. A maximum M 2 amplitude of~4 m was observed at 15 ka BP, which is likely due to the funnelling of the tidal wave through the evolving narrow channel between mainland Scotland and the Western Isles.
The degenerate amphidromic point to the southeast of Ireland was simulated to move farther north during deglaciation, as the Irish Sea flooded, and shifted towards the region of increased tidal dissipation which moved progressively eastwards after 12 ka BP to its present-day position near the eastern Irish coast (it remains degenerate). This shift in the amphidromic point will have influenced the tidal response to rising sea levels along the east coast of Ireland, which can be seen in the decreasing tidal elevation amplitudes along the eastern Irish coast since the opening of the North Channel (16e15 ka BP). These results for the northern Celtic Sea and Irish Sea differ with respect to the pattern and timings of the results presented by Uehara et al. (2006) due to the timings of submergence of the shelf in this area. In contrast to the simulations of Uehara et al. (2006), in the ROMS þ Bradley simulations, the northern Celtic Sea remained flooded from 21 ka BP to the present day with no Late-glacial land bridge between Britain and Ireland. As the semi-enclosed Irish Sea was flooded, the M 2 amplitudes remained <1 m until 16 ka BP. Between 16 ka BP and 15 ka BP, the North Channel opened, connecting the Irish Sea with the Malin Shelf/North Atlantic, and the maximum M 2 amplitude in the eastern Irish Sea doubled from~1.3 m to >2.5 m as the tide was able to propagate fully into (and through) the Irish Sea.
Along the palaeo-coastline joining Britain with northern France, M 2 amplitude reached a maximum of~4 m at 20 ka BP, and remained higher than 2.5 m until the English Channel fully connected with the North Sea between 9 and 8 ka BP. The division of this region of high M 2 amplitude into the Bristol Channel and along the Brittany coast occurred earlier in the ROMS þ Bradley simulations (after 14 ka BP) than was simulated by Uehara et al. (2006), i.e., after 12 ka BP. During this period, the areas of high M 2 amplitudes migrated along the northern coast of France as it was submerged, into the mouth of the English Channel, reaching amplitudes in excess of 4.5 m in the Gulf of St Malo (northern France) at 10 ka BP. The opening of the English Channel allowed the tide to propagate into the North Sea, reducing the tidal resonance in this region of the shelf and increasing the tidal dissipation due to friction effects in the relatively shallow channel, thus reducing the tidal amplitude. Increasing M 2 amplitudes were found in the southern Irish Sea after 14 ka BP, with the flooding of the shallow Bristol Channel. Here, the M 2 amplitude is simulated to have reached a maximum of 4.5 m at 10 ka BP, likely due to the geometry of the channel funnelling the tides, as well as due to resonant effects. In the Bristol Channel, the M 2 amplitude decreased to~4 m between 6 and 5 ka BP after which the M 2 amplitude increased to the present-day value of >4.4 m.
The modern amphidromic point in the eastern North Sea was simulated to have been degenerate between 21 and 15 ka BP (Fig. 4), with the exception of at 20 ka BP. After 15 ka BP, this amphidromic point moved progressively southeast as the North Sea deepened -at the same time the tidal amplitude along the eastern coast of Great Britain increased -and has remained in its approximate present-day location since 8 ka BP. The present-day amphidromic point in the southern North Sea/northern English Channel appeared at 8 ka BP, which is 1 ka earlier than was simulated by Uehara et al. (2006), and shifted eastwards to its presentday position by ca. 6 ka BP. As found by Uehara et al. (2006), there were no significant changes in M 2 amplitude from 6 ka BP to present.

Tidal-induced bed shear stress
Large-scale redistribution of sediments by hydrodynamical processes occurs in shelf seas, influencing basin and coastal evolution. Shelf sea sediments are important, as they comprise the sea bed, determine the turbidity of water, provide a substrate for marine benthic organisms, carry organic matter, and are involved in biogeochemical exchanges. In the majority of shelf sea and coastal regions, both waves and currents play a role in sediment dynamics; however, the combined effect is not simply a linear addition of the two independent effects (Soulsby, 1997). The frictional force exerted by currents on the seabed is referred to as bed shear stress (t 0 ) and is expressed as the force exerted by the flow per unit area of bed in terms of the density of water (r) and the frictional velocity (or shear velocity, u * ), i.e., t 0 ¼ ru 2 Ã . Sediment transport (of noncohesive sediments) occurs when the bed shear stress exceeds the threshold of motion, which is a dimensionless form of the bed shear stress; hence, in several modelling studies of the region, the directions of peak bed shear stresses have been used as a proxy for net sediment transport (e.g., Pingree and Griffiths, 1979;Aldridge, 1997;Hall and Davies, 2004;Neill et al., 2010;Ward et al., 2015).
At present, regions of particularly high bed shear stresses are found in the Irish Sea and the North Channel, along the English Channel, near the Wash (east coast of Great Britain) and around the Orkney Isles (Scotland) and Faroe Isles (Denmark). Such regions of high bed shear stresses are often coincident with coarser sea bed sediments (Holmes and Tappin, 2005;Ward et al., 2015), as is reflected in the British Geological Survey seabed sediment map 1 . The regions of simulated high peak bed shear stress on the northwest European shelf seas have migrated considerably since 21 ka BP (Fig. 5), resulting from the increasing water depths and the changes in areal extent of the shelf seas. The focus of the results is on the bed shear stresses output by the new ROMS þ Bradley simulations; however, the differences in simulated bed shear stresses between the three palaeotidal models, for the 12 ka BP time slice, are illustrated in Fig. 6. This figure shows that, despite the different hydrodynamic models used in KUTM þ Lambeck (2D) and ROMS þ Lambeck (3D), and in the way the bed shear stress is calculated (2D vs 3D), the simulated bed shear stresses are much closer in value than the ROMS þ Bradley bed shear stress values.
Between 21 and 12 ka BP, there was a broad region of high peak bed shear stresses in the southern Celtic Sea. At 21 ka BP, the water depth was shallower than 70 m, and a maximum bed shear stress of 16 N m À2 was simulated, which relates to a current speed of 2.3 m s À1 through the peak bed stress t ¼ rC D ujuj, where C D is the drag coefficient (0.003), r is the density of seawater (1023 kg/m 3 ), and u is the velocity vector. This region of high bed shear stresses in the southern Celtic Sea migrated northeast (onto the shelf) with rising eustatic sea levels, into the English and Bristol Channels as they formed. A clear divide in high bed shear stress occurred in this region after 10 ka BP, moving into the English Channel and the Irish Sea. With these palaeobathymetries, between 16 and 15 ka BP, the North Channel opened as the gateway between the Irish Sea and the North Atlantic . This resulted in higher simulated bed shear stresses in the Irish Sea after this time, peaking at~18 N m À2 (~2.4 m s À1 bottom current speed) at 12 ka BP. Before 15 ka BP, high bed shear stresses in this region were not observed in the semi-enclosed Irish Sea, which was characterised by fairly low tidal current speeds (Fig. 5). By comparison, the Lambeck GIA model predicts the North Channel opened later, between 15 and 14 ka BP.
Very high peak bed shear stresses were also simulated at 21 ka BP south and west of the Faroe Isles (maximum of 24 N m À2 ), and remained high until 16 ka BP. The minimum present-day water depth for the same region south of the Faroes is 120 m (in comparison to <20 m at 21 ka BP), and the peak bed shear stress in the region is an order of magnitude lower (maximum of 2.4 N m À2 ). There was also an area of bed shear stress divergence northwest of the Faroe Isles, most apparent before 14 ka BP, across the Faroe-Scotland ridge. The area of high bed shear stresses north of Scotland (around what is now the Orkney archipelago) has decreased since the region was first submerged between 16 and 15 ka BP. There are still particularly high bed shear stresses south of Orkney (the Pentland Firth), where peak simulated present-day tidal current speeds are >3.5 m s À1 (i.e., bed shear stress >34 N m À2 ).
Between 21 and 12 ka BP, there were high bed shear stresses along the western shelf break, in particular along the edge of the Malin Shelf and west of Ireland, due to the flooding of the shelf in this region. The Bradley palaeotopography for 21 ka BP suggests the British-Irish Ice Sheet extended out to the shelf break in this area. As the North Sea flooded, an area of high bed shear stresses emerged after 14 ka BP, southwest of Dogger Bank and along the north of the palaeo-coastline. Bed shear stresses have remained high along the east coast of the Great Britain at this latitude, in the area of the Wash.

Constraining sea level index points (SLIPs)
Models of GIA are continuously being revised and improved as more observational data of increasing quality become available for constraining the models, and also as understanding of the GIA and glaciation-deglaciation processes improves (Peltier, 1994;Lambeck, 1996;Lambeck et al., 1998;Milne et al., 2006Milne et al., , 2009Shennan et al., 2000b;Peltier, 2002a;Brooks et al., 2008;Bradley et al., 2009Bradley et al., , 2011Teferle et al., 2009). Such observational data, referred to as sea level index points (SLIPs) provide the primary observational evidence for constraining GIA models. Knowledge of the age, location, altitude and tendency (i.e. the level of marine influence) of such sediments or features allows the records to be used as SLIPs (Shennan et al., 2006a). Because few SLIPs formed exactly at palaeo-mean sea level, the 'indicative meaning' of a SLIP relates its vertical relationship to a former tide level using two parameters (Shennan et al., 2006a). The first of the two parameters is the reference water level (e.g. mean high water of spring tides, MSL þ M 2 þ S 2 tidal amplitude), which allows for calculation of the sea-level change relative to present by comparing the reference water level and the altitude of the reference point to the national datum. The second parameter considered is the indicative range, which is the vertical range over which the accumulating sediment could occur (Shennan et al., 2006a). MHWST varies with tidal amplitude, which is influenced by RSL changes (Gehrels et al., 1995;Shennan et al., 2000a;Uehara et al., 2006;Neill et al., 2010), therefore correcting SLIPs for palaeo-changes in MHWST is critical for constraining their indicative meaning. Fig. 7 shows the changes in simulated MHWST at Arisaig, on the west coast of Scotland, which is the location of the most extensive SLIP record for the British Isles (Shennan et al., 2006b), with data covering the last 16 ka. This figure revises the MHWST evolution curve presented by Neill et al. (2010), which was developed using the GIA model of Lambeck (1995), and the POLCOMS tidal model (3D, 1/12 ). Between 18 ka BP and 8 ka BP, there were considerable differences between simulated tidal amplitudes and RSL. During early deglaciation, the large tidal amplitudes were due to the large tidal range along the western extent of the British-Irish Ice Sheet (Fig. 4). The rapid decrease in RSL (i.e., sea level fall) during early deglaciation was due to the rapid GIA caused by the decrease in loading from the diminishing British-Irish Ice Sheet.

An updated palaeotidal model for the northwest European shelf seas
The numerical modelling of both tides and GIA has progressed significantly in recent years with advances in observational and computational techniques, and as such there is the need for new palaeotidal simulations to reflect these recent developments. A new, high resolution (~1/24 ), three-dimensional (3D) palaeotidal model for the northwest European shelf seas has been developed (ROMS þ Bradley). The new palaeotidal model incorporates the latest information on dynamic palaeotopography for the region . Prior to this work, the most recent palaeotidal models of the region include those by Uehara et al. (2006) (2D, ROMS þ Bradley) and Neill et al. (2010) (3D), both of which were 1/ 12 (~10 km) spatial resolution, and so did not fully resolve many features of the shelf seas. A second new palaeotidal model has been developed (ROMS þ Lambeck), which incorporated the older GIA model of Lambeck (1995). In the Bradley and Lambeck GIA models, the global ice sheet model was relatively low resolution, at 1/2 in longitude and 1/4 in latitude, respectively. The 25 km spatial resolution ice model GB-3 developed by Lambeck (1993) was used as the regional ice model in the Lambeck model. The local ice sheet model used by Bradley et al. (2011) was constrained using reconstructions by Shennan et al. (2006b) and Brooks et al. (2008). For the non-local ice-melt signal, the global ice model of Bassett et al. (2005) was used, with the British-Irish Ice Sheet component removed (to allow for the new local ice sheet model of Bradley et al. (2011)). The Bradley model takes into account considerably more field data regarding the extent of the ice sheet which were not available when the Lambeck model was developed. Considerable differences were found between the palaeotidal models of the region which incorporated the GIA models of Lambeck (1995) and Bradley et al. (2011) (see Table 1 for model details), although not between different hydrodynamic models whose palaeotopographies were derived from the same GIA models (i.e. ROMS þ Lambeck, KUTM þ Lambeck). The inter-model comparison presented here indicates that the main influence on the timings and magnitude of changes to tidal dynamics is RSL (Fig. 6).
Another important consideration when making inter-model comparisons is whether model outputs are 2D or 3D, and hence whether tidal currents are calculated using bottom or depth- averaged tidal current speeds, which has implications for simulated bed shear stresses. For example, for considering tidal-induced sediment transport, the modelled bottom bed shear stress is favoured over bed shear stress calculated using depth-averaged tidal currents speeds, since the latter make more significant assumptions about the vertical distribution of velocity. The advantage of the palaeotidal models developed using ROMS is the versatility offered by variables being output in both 2D and 3D. With these points in mind, and considering the new palaeotidal model developed here has more extensive output than that presented by Uehara et al. (2006), the ROMS þ Bradley model is considered to be a significant improvement on existing palaeotidal models of the northwest European shelf seas.

Regional implications of the new palaeotidal model
The direction of peak bed shear stress can be an indicator of sediment movement (Pingree and Griffiths, 1979;Aldridge, 1997) and the magnitude and direction of bed shear stresses have evolved with changing sea levels (Hall and Davies, 2004;Uehara et al., 2006;Neill et al., 2010). Understanding long-term and large-scale sediment dynamics and geological processes that have formed coasts and the surrounding bathymetry is important for the appropriate management of the present-day coastline as well as for predicting the impact of future climate and sea level change. For example, the tidal elevation amplitudes output by ROMS þ Bradley in the northwest Celtic Sea, in the vicinity of linear tidal sand ridges (e.g., Belderson et al., 1986;Scourse et al., 2009b), is considerably smaller during early deglaciation than those modelled by Uehara  (2006). The northern extent of the sand ridges extend into this area of lower elevation amplitudes, which had correspondingly lower bed shear stress. The simulations from ROMS þ Bradley thus pose a problem regarding defining the origin of the sand ridges, which are hypothesised by Scourse et al. (2009b) to have formed by high tidal-induced bed shear stress remobilising sediments between 20 and 12 ka BP (or until 10 ka BP in the southerly region). An alternative hypothesis was proposed by Praeg et al. (2015), who suggested that the sand ridges in this region are not in fact tidallyreworked sediment features, but are relic features from the advance and retreat of an Irish Sea Ice Stream which extended farther south than previously proposed limits (e.g. Scourse et al., 1990;Svendsen, 2004;Sejrup et al., 2005). Outputs from KUTM þ Lambeck have previously been used to support empirical studies, such as the evolution of the Celtic Sea mega 'linear tidal sand ridges' (Scourse et al., 2009b) and the formation of sediment waves in the Irish Sea (Van Landeghem et al., 2009a). We have shown here that the predicted bed shear stress direction is sensitive to the GIA model used, the same finding was reported by Uehara et al. (2006). We find here that with the ROMS þ Bradley simulations, the direction of the peak bed shear stress vector in the southern Irish Sea was found to reverse between 11 and 10 ka BP. This is earlier than was described by Van Landeghem et al. (2009a), who found that use of the less realistic 'fixed' boundary conditions of the KUTM þ Lambeck model actually provided a better fit with observed directions of bottom bed shear stress vectors in the Irish Sea (15 ka BP to present day). The bed shear stress vectors in the southern Irish Sea were directed southwest for a longer time period in the new palaeotidal model (ROMS þ Bradley), which is more consistent with the observed bedload transport direction, since Van Landeghem et al. (2009a) reported opposing observed and predicted bed shear stress vector directions in this region before 9 ka BP.
In addition to these differences in the bed shear stress magnitudes, the area largely remained marine throughout deglaciation in the ROMS þ Bradley simulations, whereas the RSL was considerably lower in the Celtic Sea in the ROMS þ Lambeck and KUTM þ Lambeck simulations. From observations, the grounding line of the Irish Sea Ice Stream is thought to be farther south into the Celtic Sea (Scourse et al., 1990;Chiverrell et al., 2013;Praeg et al., 2015) than is represented by these palaeotopographies. This is likely due to the lowresolution ice sheet models used here (combination of Bassett et al., 2005;Shennan et al., 2006b;Brooks et al., 2008) for defining the location of the ice-water boundary in the palaeobathymetric grids in each of the palaeotidal models.

Constraining sea level index points (SLIPs)
Constraining SLIPs with MHWST corrections is a somewhat iterative process, in which SLIPs are used to constrain a GIA model, the output of which is used as input for a palaeotidal model, which itself is necessary for correcting the SLIPs from their in-situ position. As proposed by Neill et al. (2010), one approach to this problem would be using split calibration verification methods, whereby only a portion of the SLIPs are used to constrain GIA model output, which are used to construct palaeobathymetric grids for the palaeotidal model (this is standard practice for many model calibration/validation exercises). A tidal correction could then be applied to the remaining portion of the SLIPs using the palaeotidal model output, and the corrected SLIPs could be used to further constrain the GIA output. Regardless of this circularity, the outputs of this palaeotidal model are better-suited for such application, since they are higher resolution than previous palaeotidal models (and so the locations of SLIPs can be defined with more accuracy), and incorporate an improved GIA model for the British Isles. For example, Fig. 7 highlights the differences between the modelled MHWST (tidal amplitudes) at the nearest grid cell to Arisaig for the different palaeotidal models, and which is the correction that could be applied to the Arisaig SLIP to adjust for evolving changes in tidal amplitudes. The actual distance from the location of Arisaig varies with model resolution, with the modelled land/ice and water mask, and with different GIA models. Again, it is assumed that ROMS þ Bradley is the preferred palaeotidal model, since it is of higher spatial resolution (1/24 ), and incorporates the latest advancements in GIA modelling for the region. This revised sea-level curve updates the MHWST evolution curve presented by Neill et al. (2010) and brings the corrected SLIP observations into closer agreement with the GIA predictions for this area.

Future work
Although the new ROMS þ Bradley palaeotidal model is a significant improvement over previous palaeotidal models of the region (e.g. Uehara et al., 2006;Neill et al., 2010), there is potential to further develop the simulations by increasing the temporal and spatial resolution, depending on the availability of boundary tidal forcing. For example, the grid resolution could be further refined, and the simulations run for 500 year time slices, to enhance the simulation of evolving tidal dynamics in specific regions, or during particular periods of interest, such as when there are large temporal gradients. Alternate wetting and drying should be included in simulations with increased spatial resolution, as there would be increased need to resolve the intertidal regions. For more comprehensive palaeotidal modelling, future models should be forced with more tidal constituents than those considered within this study, where available. Any future palaeotidal models of the region should include any recent developments in GIA modelling (e.g. Peltier et al., 2015), which themselves incorporate improvements in sea-level models and more accurate ice sheet reconstructions and new data on RSLs (e.g. Chiverrell et al., 2013;Sturt et al., 2013;Praeg et al., 2015).
Only tidal-induced currents and bed shear stresses are considered within this paper. Although long-term sediment dynamics on the northwest European shelf seas are dominated by tidal currents, waves also contribute significantly (e.g. van der Molen, 2002;Neill et al., 2009). Waves are the primary mechanism for inter-annual variability, due to sensitivity to variability in atmospheric (wind) forcing (Neill et al., 2009). Thus, a significant, yet challenging development in this field of palaeo-modelling would be to conduct coupled tide-and wave modelling (e.g. Warner et al., 2010). Coupled tide-and wave modelling can be computationally expensive, but it would increase the accuracy of the simulations of changing hydrodynamics, in particular where bed shear stresses are of interest, since these coupled model outputs would include consideration of wave-induced currents. A significant limitation of this coupled model approach is the lack of suitable palaeowind data for the northwest European shelf seas, which is required as boundary forcing for palaeowave models. For example, a previous palaeowave modelling study of this region (12 ka BP -present), conducted by Neill et al. (2009), used a range of present-day wind fields as model forcing as no suitable palaeowind field were available. Future work should focus on the development of higher resolution palaeowind fields for the region, at suitable time slices for progressing coupled palaeotide-and wave modelling.
Furthermore, finding a robust method for constraining and validating palaeotidal models across domains using observational data remains to be addressed, as the only existing data point for validating palaeotidal models was developed by Austin and Scourse (1997) and Scourse et al. (2002), was used by Uehara et al. (2006) to compare the skill of different palaeotidal models. This data point was based on the timing of stratification onset in the Celtic Deep, and was limited to a single sediment core location.

Conclusions
A new 3D palaeotidal model (ROMS þ Bradley) for the northwest European shelf seas has been developed which incorporates the most up-to-date GIA model for the region. By comparing the outputs of simulated tidal dynamics with those from other palaeotidal models for the region (ROMS þ Lambeck and KUTM þ Lambeck), which incorporated a less well-constrained GIA model, we show that of most influence on the model outputs is the GIA model choice, rather than the hydrodynamic model choice (2D versus 3D model). Although the versatility (and comprehensiveness) of a model which outputs in both 2D and 3D is preferable over a depthaveraged-only tidal model, this is not the key consideration with regards developing palaeotidal models. The differences in the tidal dynamics simulated using different palaeotopographies is also significantly more than the variability induced by different hydrodynamic models. Thus, use of the most up-to-date and well constrained GIA model is important in order to reproduce past RSLs as accurately as possible. With this in mind, and considering the palaeotidal model developed here has more extensive output than that presented by Uehara et al. (2006), this new palaeotidal model (ROMS þ Bradley) is an improvement on existing palaeotidal models. Work on constraining GIA models to produce better fit with observational data on RSLs is ongoing, as significant uncertainties remain in the inputs (i.e., the ice-, earth-and sea-level models) to GIA models.