Early Holocene Greenland-ice mass loss likely triggered earthquakes and tsunami

a Lantmäteriet, Lantmäterigatan 2, 80182 Gävle, Sweden b Department of Geosciences, Virginia Tech, 4044 Derring Hall, Blacksburg, VA 24061, USA c Center for Coastal Studies, Virginia Tech, 926 West Campus Drive, Blacksburg, VA 24061, USA d Department of Physics and Physical Oceanography, Memorial University of Newfoundland, St. John’s, Newfoundland A1B 3X7, Canada e Department of Earth and Environmental Sciences, University of Ottawa, STEM Complex, 150 Louis Pasteur Private, Ottawa, Ontario K1N 6N5, Canada f Department of Geography, Durham University, Lower Mountjoy, South Road, Durham, DH1 3LE, UK g Geological Survey of Denmark and Greenland, Øster Voldgade 10, 1350 Copenhagen K, Denmark


Introduction
The growth and decay of ice sheets during a glacial/interglacial cycle affect a multitude of processes on the surface as well as in the interior of the Earth, which are commonly termed glacial isostatic adjustment (GIA). For example, the mass redistribution of water between the ice sheets and the oceans causes changes in the Earth's shape, gravity, rotation and sea level (Wu and Peltier, 1982). This climate-driven surface loading results also in significant horizontal and vertical stress changes due to the enormous mass of the ice sheets (Johnston, 1987). In a compressional stress setting, where horizontal stresses are larger than the vertical stress, fault slip and thus earthquake activity is inhibited (Johnston, 1987), hence partially explaining the relatively low seismic activity in present-day Greenland (Voss et al., 2007). During deglaciation, however, the vertical stress decreases in relation to the vanishing The ice sheet undergoes negative mass balance in response to climate warming. (B) Ice sheet retreat causes a viscoelastic glacial isostatic response from the solid Earth. (C) Due to an asynchronous decrease of horizontal and vertical stresses in a compressional stress setting, a pre-existing fault is reactivated triggering an earthquake and tsunami. Khan et al., 2016). We present the first calculations of stress changes induced by ice-mass loss in the early Holocene to assess whether and where glacially-triggered earthquakes were likely to have occurred. Our modelling results indicate that southern Greenland is the most prone to glacially-triggered faulting. Using a faulting scenario that is consistent with geological records of relative sea level change, we show that a modelled fault could have been reactivated offshore, thus producing a tsunami wave (Fig. 1C) during the early Holocene for which we compute the wave height distribution around North Atlantic coasts. As this analysis consists of a combination of several different methods, each producing results required by another, to avoid repetition we include a description of the methods within an extended Section 2 on Methods, Results and Discussion.

Stress modelling
A recent model reconstruction of the Greenland Ice Sheet (termed Huy3; Lecavalier et al., 2014) and its appendant (optimal) Earth viscosity model (lithospheric thickness of 120 km, upper mantle viscosity of 5 × 10 20 Pas, lower mantle viscosity of 2 × 10 21 Pas) and ocean-load model are used to calculate stress changes during the past 120,000 years for the whole of Greenland in a compressional stress setting. Our modelling procedure follows the approach described in Wu (2004) which uses a three-dimensional (3D) flat Earth model using the finite-element software ABAQUS (Hibbitt et al., 2016) to estimate GIA-induced displacements and stresses. The model consists of eleven layers with different material parameters covering a depth range from the Earth's surface to the core-mantle boundary. The upper four layers are purely elastic and form the lithosphere, while the lower seven layers represent the visco-elastic mantle. The horizontal length scale of the finite elements is 50 km and the vertical length scale gradually increases from 5 km in the crust (upper 30 km of the lithosphere) to a maximum of 590 km in the lower mantle. No lateral variations of the material parameters are used. Earth model variations have only a small effect on the reactivation time for faults located between the ice margin and the ice-sheet centre (e.g., Kaufmann et al., 2005), and the usage of a 3D Earth model is therefore not needed here. The model domain features Greenland in its centre and is square in shape with a side length of 4500 km. The model domain also includes parts of North America and northern Europe which are, as in Lecavalier et al. (2014), loaded with a North American ice sheet model by Tarasov et al. (2012) and a Fennoscandian and Iceland ice sheet model by Peltier (2004) to incorporate the GIA response from these regions as well. The Huy3 ice model has a horizontal resolution of 20 by 20 km in Greenland, which has been upscaled to 50 by 50 km for the implementation on the finiteelement grid. Ocean mass changes are not calculated within the GIA model directly as the here applied flat model approximation is not able to solve the sea-level equation. Thus, an ocean load obtained from a 1D spherical GIA calculation (Mitrovica et al., 1994) is used, which is based on the same ice and Earth model configuration as the finite-element model. The ice and ocean load are then applied together to the Earth model to obtain GIA-induced displacements and a stress field. However, the obtained stress field is not complete as a different equation of motion is used in the finiteelement methodology in ABAQUS (Hibbitt et al., 2016) as required for GIA purposes. Thus, the obtained stress tensor is transformed to GIA-induced stresses (Steffen et al., 2015) while the displacement vector is the same. The calculated GIA stresses are then combined with tectonic background stresses to analyse the potential for GIF reactivation (Steffen et al., 2014).
Compressional stresses are applied to simulate a ridge push force east of Greenland at the Atlantic mid-ocean ridge acting towards the stable craton of North America to the west (Bird, 2003) and resulting in maximum east-west horizontal stress and minimum north-south horizontal stress. This background stress field for Greenland is confirmed by mantle flow models (Conrad and Behn, 2010), which indicate a horizontal mantle movement for the area offshore southern Greenland. In addition, a horizontal direction of the maximum stress is inferred in the World Stress Map offshore of the north-east of the United States (Heidbach et al., 2018), and studies of the palaeostress direction in the Palaeocene show a rotation of the maximum horizontal stresses from north-south to east-west going from the north-eastern coast (Peary Land) to the south-eastern coast (Skjoldungen; Guarnieri, 2015). Stresses in the north-south direction are likely to be small and similar in magnitude to vertical stresses. This is due to the lack of active plate boundaries to the south and north of Greenland (Bird, 2003), and is further supported by observations of palaeostresses along the eastern coast, which show a decrease of the intermediate stress magnitude from north to south (Guarnieri, 2015). Such a decrease in the intermediate stress magnitude relates to an increase of the stress ratio R from north to south (Guarnieri, 2015). This parameter links the maximum, medium and minimum stresses with each other (Etchecopar et al., 1981), and a high stress ratio of 0.95 is used here.
The background stresses are calculated separately and not modelled as part of the GIA model, therefore no plate boundaries need to be explicitly included into the GIA model. The principal stresses of the background stress field are determined based on the following equations: with ρ as the density (constant density is assumed), g as the gravitational acceleration, h as the depth, and μ as the coefficient of friction (0.6 is assumed here). The equations are based on the relation between maximum and minimum principal stress defined by Anderson (1951) and are valid for a thrust-faulting background stress regime. The principal stress directions are aligned with the x-, yand z-directions (σ 1 x east-west, σ 2 y north-south, σ 3 z depth). The estimated stress magnitudes are then combined with the glacially-induced stress field to find the change in the Coulomb Failure Stress.

Coulomb failure stress changes
We use critically stressed conditions in the crust, which is valid for intraplate areas (Zoback and Townend, 2001), and analyse the change in Coulomb Failure Stress CFS (Harris, 1998) which helps visualize stable and unstable seismic conditions. Put simply, positive values of this quantity represent unstable conditions indicating that seismic activity is likely, and negative values point to stable conditions where seismic activity is unlikely. Only two areas, the southern tip and northern coast of Greenland, experience unstable conditions in the early Holocene due to ice retreat (Fig. 2, Movie S1).
Seismic activity within deglaciating regions requires pre-existing faults, which can be reactivated to release the deglaciation-related stress build-up (Steffen et al., 2014). Faults in North Greenland strike mainly east-west (90 • /180 • ) while those in southwestern Greenland, close to the small town of Nanortalik, strike mainly northwest-southeast (135 • /315 • ) and northeast-southwest However, while glacially-triggered earthquakes in northern Europe have been identified using topographical changes and visible fault outcrops (Lagerbäck and Sundh, 2008), such large-scale features have not been observed in southern Greenland. Therefore, we look to observations of past relative sea-level (RSL) change to determine whether they are compatible with the timing and amplitude of faulting suggested by our model.

Relative sea level history in Nanortalik, Southern Greenland
Several RSL data sets have been collected in the area around Nanortalik ( Fig. 3) with sea-level indicators from 13.511 ± 0.236 ka BP onwards (Bennike et al., 2002), thus covering the time when the area is predicted to have become unstable (Fig. 2D). These RSL data are of the highest quality -they are based on 14 C dated sediment samples from isolation basins which are rock depressions in the landscape that have been uplifted and isolated from the sea in the past (Bennike et al., 2002). The Nanortalik RSL data show rapid early Holocene sea-level fall from at least ∼32 m above present around 13.8 ka cal BP, reaching present-day RSL by c. 10 ka cal BP and continuing to a lowstand in the early Holocene before rising to present in the late Holocene (Fig. 3C). Other Holocene RSL data in this region from Qaqortoq, 90 km NW of Nanortalik (Sparrenbom et al., 2006), and Igaliku, 100 km N of Nanortalik (Bierman et al., 2018) also show rapid early Holocene RSL fall.
We apply the deglaciation history of Huy3, which is based on a Greenland-wide ice extent and RSL database , alongside the accompanying 1D Earth model within a GIA model to investigate the fit of the Huy3 GIA model to the Nanortalik RSL data. The majority of RSL data from around Greenland can be explained by Huy3 model reconstructions . However, a misfit exists for data points in southern Greenland, and particularly at Nanortalik, especially for the four oldest data points. The Huy3 model resulted from extensive sensitivity tests that explored various atmospheric and oceanic (sea-level) forcings, and ice sheet model parameters to best capture the observational constraints. This sensitivity analysis also included an exhaustive sampling of various spherically symmetric Earth model parameters (see Fig. 10 in Lecavalier et al., 2014): lithospheric thicknesses ranging from 71 to 120 km, upper mantle viscosities ranging between 1 × 10 20 to 5 × 10 21 Pas, and lower mantle viscosities spanning 1 × 10 21 to 50 × 10 21 Pas. This extensive parameter search was unable to resolve the first-order misfit to the Nanortalik data ( Fig. 3; Lecavalier et al., 2014). A more recent study, using a 3D Earth model in combination with the Huy3 ice model was not able to improve the data-model misfit at Nanortalik by considering a number of different estimates of lateral viscosity variations (Milne et al., 2018). While the addition of lateral structure improved the data-model fits in the early Holocene, those for later times were made worse (Fig. 3C). Two additional ice model reconstructions were considered (ICE-5G (Peltier, 2004) and ANU (Fleming and Lambeck, 2004)) but the model fits were of lower quality than those for the Huy3 model (Fig. S2) and even larger discrepancies exist. We therefore excluded these models for the estimation of the glacially-induced stress field.
RSL predictions from Huy3 with the new 3D Earth model indicate the deglacial marine limit was reached at c. 11 ka cal BP and rapid RSL fall occurred immediately thereafter (Fig. 3C). The ice history here is most likely inaccurate as suggested by  and Milne et al. (2018) because the timing of the marine limit being reached, which should correspond to initial deglaciation at the location, is ∼3 ka too late given the evidence of ice-free conditions at lake N14 at 13.8 ka cal BP. Despite this issue with the timing of deglaciation, and therefore the timing of initial RSL fall in the Huy3 predictions, there still remains a significant discrepancy between the elevation of the RSL data before ∼10.6 ka cal BP and the GIA model predictions. In particular lakes N14 and N18 are up to 14 m above the upper limit of uncertainty in the 3D Earth model (and 24 m above the lower uncertainty in this model) (Fig. 3C), and lakes N19 and N24, whilst falling within the upper limit of the 3D Earth model uncertainty, currently suggest a faster initial RSL fall compared to what might be predicted by an ice model with earlier deglaciation in this region .
We therefore propose the hypothesis that tectonic activity led to the movement of the four RSL index points at Nanortalik older than ∼10.6 ka. The occurrence of such an event would influence the elevation of RSL index points older than this age but not younger ones, thus bringing the RSL data into closer alignment with the Huy3 predictions (Figs. 3C, 4). As the sediments in the isolation basins N14, N18, N19 and N24 show no evidence for sea-level rise (i.e. a later transgression into the basin following initial isolation), this means a maximum correction of up to 16.5 m at any of the Nanortalik isolation basin locations is permitted by the data (Fig. 3C, red boxes). In the following section, we make a preliminary test of our hypothesis by simulating plausible faulting scenarios in southern Greenland. Invoking a faulting event dramatically improves the fit between model RSL predictions and faulting-corrected, observed RSL heights (Fig. 4).

Fault modelling
To simulate an earthquake due to the modelled stress changes, we created a two-dimensional (2D) GIA-fault model (Steffen et al., 2014) and considered a range of plausible fault parameters (e.g. fault depth, fault width, friction). The stress and displacement results obtained from the 3D GIA model and stress analysis (described above) are used on a 2D profile together with the same Earth model parameters. The 2D profile is perpendicular to the strike direction (Fig. S3) and crosses the points N14, N18 and N19, but has a distance of about 15 km to N24. The spatial resolution is greatly increased compared to the 3D model and varies between 500 m in the crust to a few kilometers in the lower part of the lithosphere. We do not apply the Huy3 ice model but use the corresponding stresses for each time step of the 3D GIA model and implement them in the 2D GIA-fault model. As dis-cussed above, the tectonic regime in southern Greenland indicates that deglaciation would most likely result in thrust faulting with a strike orientation that is NW-SE (or SE-NW). Of all the faults we considered, via a set of parameter values for dip, strike and friction (see Fig. S1), the one described below is the most likely to have been reactivated based on the offsets indicated by the RSL data. The modelled offshore thrust fault southwest of Nanortalik (Fig. 4A) results in a surface deformation that uplifts all RSL data points older than 10.6 ka in the area of Nanortalik. This fault has a dip of 45 • and extends from a depth of 5 km to 24 km and hence does not outcrop at the surface; its strike is parallel to the outer coast at 315 • and is thus parallel to surface faults identified in this area.
An additional parameter in the modelling of the fault displacement is the coefficient of friction. Steady-state and static friction values are defined along the fault surface. Observations show that The green line is obtained from the best-fitting Earth model with a lighter green bounding envelop based on RSL curves using the nominal 95% confidence interval Earth models (described in Table 2 of Lecavalier et al. (2014)). A range of RSL curve predictions using various 3D Earth models (yellow area) indicate the plausible influence of lateral Earth structure (after Milne et al., 2018). The RSL reconstructions are shown with 2-σ uncertainty in the time and height range as black bars. The maximum possible offset to allow a better fit to the RSL predictions and satisfy the observational constraint of RSL fall are shown as red squares for N14, N18, N19 and N24. the steady-state friction is about 10-30% of the static friction (Di Toro et al., 2011). The steady-state (μ ss ) and static friction (μ k ) are related to each other in ABAQUS (Hibbitt et al., 2016).
The decay coefficient d c is taken to be 1.0, as this would equal the equation presented in Di Toro et al. (2011) for laboratory earthquakes Here, we use a static friction of 0.6 and a steady-state friction of 0.12 (20% of the static friction). As the friction is unknown, we calculate an uncertainty due to this by applying steady-state frictions of 0.06 and 0.18 as well due to the observed changes between 10% and 30% (Di Toro et al., 2011), respectively. The duration of the earthquake is chosen to be 10 seconds as the average rupture velocity is between 2.6 and 3.0 km/s (Heaton, 1990), which results in a duration of 10.4 to 9 s for this specific fault with a width of ∼27 km. However, in the estimation of the uncertainty (Fig. 4A), we also consider different durations (1 s to 60 s) of the rupture propagation. Although a post-seismic phase was not included in the fault modelling, we estimated the post-seismic displacement following the earthquake using VISCO1D (Pollitz, 1997). The stress release of the earthquake is applied to the GIA and background stresses in the subsequent time steps and no further instability occurs until today at this location. We note that our calculations neglect the influence of pore-fluid pressure as no information is available on how this parameter changes in a crustal setting beneath the ice and at the ice margin during a glacial cycle.

Earthquake magnitude estimation
In our model, the fault is reactivated at 10.615 ± 0.25 ka BP with a fault slip of 43.7 m (Fig. S3), equivalent to an earthquake with a moment magnitude of about 8.3 using a fault length of 200 km. The earthquake would be followed by a post-seismic displacement of up to 0.94 m after 1 ka using the same fault configuration and earth structure parameters. The displacement at the surface is 38.1 m (Fig. 4), which would normally relate to a calculated surface rupture length of at least 800 km using standard calculation methods (Wells and Coppersmith, 1994). However, Mattila et al. (2018) showed that GIFs appear to have a higher displacement-length ratio than the standard calculation would allow, and are usually shorter than 200 km. We therefore apply a fault length of 200 km for our Greenland faulting event based on the maximum length of GIFs found in northern Europe. This creates a moment magnitude which is larger compared to estimates for other GIFs, but those previous estimates are based on fault off- sets visible at the surface and slip within the crust was most likely larger in this instance due to the increase in displacement towards the fault centre (e.g., Kobayashi et al., 2015). In addition, the ratio of average subsurface displacement to average surface displacement is mostly larger than 1 (Wells and Coppersmith, 1994) indicating a larger subsurface displacement than is observed at the Earth's surface.
Earthquakes with such a large magnitude are rare and usually occur along active, convergent, plate boundaries rather than the intraplate setting considered here. If we consider the occurrence of several earthquakes instead of one single event, the magnitude would be decreased. Recent results by Smith et al. (2018) also showed that offsets along GIFs in northern Sweden were not created in one event but rather several smaller events. For example, the total fault slip of 47.3 m could be divided into ten events with 4.73 m slip each. This decrease in the fault slip would mean smaller surface displacements, which would allow the rupture length to be reduced to 100 km (Wells and Coppersmith, 1994). This change in the fault length and fault slip results in a decrease of the moment magnitude to 5.8 for each of the ten earthquakes, which is also more commonly observed for intraplate earthquakes (e.g., Mooney et al., 2012). However, the RSL data as well as geological maps allow no differentiation between one-large magnitude 8.3 event versus a series of smaller events over decades to centuries with lower magnitudes. Thus, neither one nor several earthquakes can be excluded, with targeted field observations (e.g., fault mapping, fault dating) being necessary to identify and con-strain the occurrence and nature of seismic activity in this area during the early Holocene. Even though the strike direction of the fault (315 • ) is consistent with faults observed in the Nanortalik region, this specific fault has not been identified, which may be due to the location of the fault being offshore and the lack of highquality, high-resolution seismic data in this region (Uenzelmann-Neben et al., 2012).
The southern tip of Greenland is an area of high seismicity compared to other parts of Greenland today (Voss et al., 2016), but recent earthquake magnitudes are mostly below 3.0 (Fig. S4), which is similar to those observed today at one of the GIFs in Fennoscandia (Pärvie fault; Lindblom et al., 2015). In addition, Voss et al. (2016) noted that many more earthquakes can be identified if the seismic station density were increased (currently only two stations, Fig. S4) and the magnitude of completeness could be decreased to below 3. A large uncertainty exists also on the location of the events and is mostly above 100 km for earthquakes recorded by these two stations (Voss et al., 2016). Therefore, previous moderate to large magnitude seismicity at the location of the hypothesised reactivated fault is a distinct possibility, particularly given the large and rapid ice sheet and sea-level changes during the early Holocene.

Relative sea level data corrections
The vertical displacement at the surface induced by this modelled earthquake would increase the elevation of RSL reconstructions older than 10.6 ka by 10.8 to 19.0 m (Fig. 4A). Here we use the displacement of the co-seismic phase only as the post-seismic displacement is less than 2.5% of the co-seismic displacement. The co-seismic displacements move the observations to within the Huy3 3D Earth model uncertainty range (Fig. 4B). Thus, invoking a faulting event dramatically improves the fit between model RSL predictions and faulting-corrected, observed RSL heights (Fig. 4). The fault-corrected RSL observations include the parametric uncertainties associated with the 2D fault model as well as the unmodelled slip dependency along the fault in the strike direction. One data point, N24, is too low using the vertical displacement correction, indicating that the fault movement is excessive (Fig. 4B). Fault slip models of previous large earthquakes show strong lateral variations along strike (e.g., Kobayashi et al., 2015), which have not been modelled here. A change in the vertical displacement based on a smaller fault slip magnitude could decrease the fault offset by up to 25% over a distance of only 10 to 15 km. The location of the RSL index point N24 is 15 km along strike relative to the other RSL data. This could lead to a decrease of the fault offset from 19.0 m to 14.25 m, moving the N24 index point up to within the vertical error of stratigraphic position and within the range of possible GIA model runs (Fig. 4B). In addition, local variations in the geology and potentially a system of faults rather than just one fault could lead to slightly different fault displacements and enhance (or reduce) the quality of fit. However, this cannot be solved using the homogeneous Earth models applied here and so is a target for future research.
Although geomorphological evidence for a reactivation event c. 10,600 years ago has not been found offshore south Greenland, the misfit between the RSL reconstructions and predictions and the timing of unstable conditions does suggest that there may have been tectonic activity at this time. Importantly, seismic data offshore south Greenland have not been analysed with this event in mind. In northern Europe, GIFs have been mainly identified by visible offsets at the surface, but in recent years soft sediment structures and high-resolution elevation data have revealed new, previously undetected GIFs (Berglund and Dahlström, 2015;Smith et al., 2018). Elevation changes due to an earthquake offshore south Greenland may also be difficult to spot as the bathymetric data for this area have a poor resolution and a gradient of less than 4.5 m/km both onshore and offshore would be obtained from the vertical fault displacement at the surface (Fig. 4A), which is difficult to identify in the landscape.
As RSL data are crucial constraints for ice model calibrations we suggest that RSL data proximal to ice sheets should be investigated for vertical displacement caused by faulting due to ice retreat, otherwise ice sheet reconstructions based on RSL observations might be biased in previously glaciated regions.

Tsunami generation
Vertical displacement of the sea floor can produce tsunami waves. We use the slip distribution of the 2D model and interpolate it to a 3D distribution towards the edges of the fault (Fig. 5B). In addition, the paleo bathymetry of the Atlantic Ocean at 10.5 ka is applied instead of the modern-day bathymetry to account for the change in sea level and the displacement of the solid earth (on-and offshore). Both estimates are taken from the 1D GIA model calculation that is also used to model the RSL change. We use GeoClaw to simulate the dynamics of the generated tsunami. GeoClaw is part of ClawPack (George, 2008) and solves the depthaveraged Shallow Water Equation using a finite volume method on adaptively refining grids (Mandli and Dawson, 2014). GeoClaw has been validated and verified using the standard benchmarks as defined by NOAA (National Oceanic and Atmospheric Administration; Synolakis, 1991) and has been applied to a variety of tsunamirelated problems (Berger et al., 2011;Arcos and LeVeque, 2015). We do not consider tides in our simulations.
We first consider results for the worst-case scenario of a single event: moment magnitude of 8.3 for a fault length of 200 km and slip of 47.3 m. For this scenario, our model predicts a sizeable tsunami that would have impacted the shorelines of North American and European continents (Fig. 5A). Greenland would have experienced the largest tsunami waves (generally exceeding 1.5 m with a maximum of 7.75 m at the southern tip). North America receives tsunami waves from 1.3 m in today's Newfoundland and Labrador to about 0.3 m in the vicinity of the northeastern corridor of today's United States of America (Fig. 5A). The tsunami waves reaching Europe are up to 0.9 m in the Northern British Islands, while they slightly exceed 1.2 m on the west coast of Ireland (Fig. 5A). Along the French and Portuguese coasts, the simulated tsunami reaches 0.7 m; while the northwest coast of the African Continent exhibits maximum tsunami wave heights exceeding 0.75 m (Fig. 5A). Note that these maxima are retrieved from the simulations in 50-m water depth and the tsunami wave elevation increases as they approach the shore. This process is known as shoaling and causes the tsunami-wave amplitude to grow by a factor between 4 and 6 (Synolakis, 1991), resulting in a maximum runup of about 5.2 m to 7.8 m and 4.8 m to 7.2 m along North American and European coasts, respectively. In comparison, considering the case when the energy is released as ten separate events with fault slips of 4.73 m each along a 100 km fault (see Fig. S5), one event would have produced run-up wave heights of up to 0.85 m along the southern Greenland coast and only a few decimetres along the Canadian and European coasts (Fig. S5).
Tsunami deposits related to the offshore Nanortalik earthquake have not been identified, which suggests the scenario of numerous smaller events is more realistic. However, this might be a premature conclusion for several potential reasons. Although the far southwest Greenland coast was ice free by this time, the interaction of any tsunami waves with permanent or seasonal sea ice in coastal areas would have decreased the tsunami impact significantly. Across Baffin Bay, the Labrador coast had grounded ice extending to the present-day shoreline (Tarasov et al., 2012;Vacchi et al., 2018), as did the coast of Iceland (Peltier, 2004). Along the western Newfoundland coast, and in many other ice-free areas of the North Atlantic, RSL was metres to tens of metres below present at the time (Fig. 5A). Thus, any tsunami deposits would now lie offshore, and it is extremely unlikely that sedimentary evidence will have been preserved through the subsequent marine transgression to present (Vacchi et al., 2018). Therefore, even if a single-event earthquake and tsunami did occur it might not be possible to find related tsunami deposits. Nevertheless, given the potential hazard associated with a major tsunami at present, due to rapid and on-going mass loss of the Greenland ice sheet, we encourage field studies aimed at detecting evidence for tsunami deposits along southern Greenland and surrounding coasts during the early Holocene.

Conclusion
It is well known that the release of glacially-induced stresses leads to the creation of earthquakes, as has been documented for parts of northern Europe from geological evidence. We propose the occurrence of glacially-triggered faulting offshore Nanortalik (south-western tip of Greenland) in the early Holocene based on stress modelling and the discrepancy between GIA model predictions and RSL data from this region. The stress release could have been associated with a single, 8.3 magnitude event or a series of moderate to strong magnitude earthquakes. If the stress release was dominated by a single event, it may have generated a large tsunami with run-up heights of several metres along eastern and western North Atlantic coasts. As ice-sheet melting on Greenland is ongoing, other areas could become unstable in the future pro- viding a potential future danger for countries bordering the North Atlantic, if offshore faults were to be reactivated again. Future field studies to test our modelling results through examining geological evidence for faulting, rapid RSL changes and tsunami evidence are encouraged.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.