# Time-variable strain and stress rates induced by Holocene glacial isostatic adjustment in continental interiors

T J Craig, E Calais, L Fleitout, L Bollinger, O. Scotti

► To cite this version:

T J Craig, E Calais, L Fleitout, L Bollinger, O. Scotti. Time-variable strain and stress rates induced by Holocene glacial isostatic adjustment in continental interiors. Tectonophysics, 2023, 854, pp.229815. 10.1016/j.tecto.2023.229815. irsn-04331492

HAL Id: irsn-04331492

https://irsn.hal.science/irsn-04331492

Submitted on 8 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Time-variable strain and stress rates induced by Holocene glacial isostatic adjustment in continental interiors

T.J. Craig<sup>a, *</sup>, E. Calais<sup>b, c</sup>, L. Fleitout<sup>b</sup>, L. Bollinger<sup>d</sup>, O. Scotti<sup>e</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

### Keywords:

Postglacial deformation

Intraplate deformation

Continental seismicity

Strain-rates

## ABSTRACT

In continental interiors, tectonically-driven deformation rates are low, often to the point where they are undetectable with modern geodesy. However, a range of non-tectonic surface processes, particularly relating to hydrological, cryospheric, and sedimentological mass changes, can produce strain-rates which on geologically-short timescales are substantially greater than those produced by tectonics. Here, we illustrate the problem that such transient strain rates may pose in low-strain environments by considering the impact that the growth and decay of the Fennoscandian and Laurentian ice sheets over the Holocene had on Europe and North America respectively. Induced deformation extended far beyond the periphery of the ice sheets, with the potential to impact on seismicity rates thousands of kilometres south of the maximum ice extent. We consider how the modelled non-tectonic deformation would have interacted with several known active fault systems, including the European Cenozoic Rift System and the New Madrid fault system. In low strain continental interiors, seismic hazard assessment – crucial for the long-term planning of critical infrastructure, including nuclear waste disposal – is often dependent on sparse information from observational and historical seismicity, and from paleoseismological studies of surface fault systems. We recommend that for a more complete seismic hazard assessment, the impact of non-tectonic transients should be considered – both in the context of the role such transients may have played in recent seismicity, and the role they may play in seismicity to come. Whilst such consideration has previously been given to the direct impact on glacial loading in areas directly glaciated, we show that it should also be considered much more broadly.

## 1. Introduction

The sparse distribution and often clustered occurrence of large earthquakes in slowly-deforming plate interiors challenges our understanding of the underlying causes of such seismicity, and hampers efforts to reliably determine the seismic hazard in these areas (e.g., Camelbeeck et al., 2007; Calais and Stein, 2009; Stein and Liu, 2009; Hough and Page, 2011; Liu and Stein, 2016; Calais et al., 2016). Modern space geodesy remains unable to detect the localised build up of elastic strain around faults in continental interiors, even in areas where large earthquakes have repeatedly occurred (e.g., Calais et al., 2005; Craig and Calais, 2014; Boyd et al., 2015). As a result, seismic hazard assessment for such areas relies on historical and instrumental seismicity catalogues

and, where available, paleoseismic studies of active fault systems. However, in such slowly-deforming regions, seismicity catalogs only capture a short-duration time interval of the fault activity, and are unlikely to be representative of their longer-term seismogenic potential (e. g., Stein et al., 2012).

In addition, the usual assumption that paleo-earthquakes, when they can be identified and characterised, occurred under strain rates that are equivalent to the present-day ones – and are therefore relevant guidelines for short-term hazard assessment – may not be valid (Craig et al., 2016). Indeed, contrary to plate boundary settings where interseismic strain rates are largely dominated by tectonic loading, strain rates in plate interiors can be significantly affected by transient non-tectonic processes that overwhelm the very slow – if any – tectonic loading.

<footnote><sup>a</sup> COMET, Institute of Geophysics and Tectonics, School of Earth and Environment, University of Leeds, Leeds LS2 9JT, UK</footnote>

<footnote><sup>b</sup> École normale supérieure, Department of Geosciences, Université PSL, 24 rue Lhomond, 75231 Paris, France</footnote>

<footnote><sup>c</sup> Université Côte d'Azur, CNRS, IRD, Observatoire de la Côte d'Azur, Géoazur, France</footnote>

<footnote><sup>d</sup> CEA, DAM, DIF, 91297 Arpajon, France</footnote>

<footnote><sup>e</sup> IRSN/PRP-DGE/SCAN/BERSSIN, 92262 Fontenay-aux-Roses, France</footnote>

<footnote>* Corresponding author.
E-mail address: t.j.craig@leeds.ac.uk (T.J. Craig).</footnote>
Examples abound of changes in surface or near-surface loading that result in measurable deformation of the lithosphere, with the potential to influence seismicity (e.g. Muir-Wood et al., 1989; Heki, 2003; Mazzotti et al., 2005; Luttrell et al., 2007; Bettinelli et al., 2007; Lagerbäck and Sundh, 2008; Calais et al., 2010; Karow and Hampel, 2010; Amos et al., 2014; Craig et al., 2016; Craig et al., 2017; Johnson et al., 2017; Rollins et al., 2020). Such load changes can result from a number of causes acting over a range of timescales, from the annual and sub-annual variation of seasonal hydrological loads, to the kyr-timescales of ice sheet variations, or to the Myr-timescales of large-scale sediment removal and redistribution. Similarly, they can operate at a variety of spatial scales, from the relatively localised deformation that results from the anthropogenic removal of groundwater, or the modulation of local surface loads caused by the volume change of major lakes, to the continental scale of major ice sheets, or the global effect of changing ocean volumes.

Whilst at plate boundaries, and in regions of relatively rapid tectonic deformation, the rates of deformation induced by such surficial processes are typically swamped by the underlying tectonically-driven deformation, in slowly deforming plate interiors the deformation rates driven by surface processes may in contrast be far greater than any underlying tectonic signal. This can result in a strain-rate field that is dominated by short-term transients, and may not, at any given point in time, be representative of the underlying stress or strain state of the crust, or of the longer-term trend in strain accumulation. A classic example is the dominant influence of post-glacial rebound in the present-day geodetic strain-rate field of tectonically-stable central-eastern North America and Fennoscandia (Nocquet et al., 2005; Calais et al., 2006; Sella et al., Jan 2007; Kierulf et al., 2014; Kreemer et al., 2014; Kreemer et al., 2018). In areas where such a non-tectonic overprint is present – or has been present over the timescales used in paleoseismological studies – one must be cautious equating strain release by paleoearthquakes to present-day strain (or stressing) rates on faults. The extreme case for this is in Fennoscandia, where the crust overlain by major icesheet thicknesses during the Last Glacial Maximum (LGM hereafter) is well-established to have hosted a number of major active faults and inferred earthquakes over the 10 ka since the last decay of the icesheet (e.g., Muir-Wood et al., 1989; Wu et al., 1999; Lagerbäck and Sundh, 2008; Craig et al., 2016; Ojala et al., 2019).

Much of continental Europe, with the exception of the Alpine orogenic belt and the Balkans, is commonly regarded as a stable continental interior, characterised by low levels of seismic activity. Geodetically observable strain accumulation related to ongoing tectonic deformation is yet to be conclusively detected (Nocquet, 2012), but is likely to be < 1 × 10⁻⁹ yr⁻¹ across the continental interior. However, major earthquakes have occurred sporadically (e.g., Basel, 1356; Dover Strait, 1580; Verviers, 1692; Düren, 1756; Lisbon, 1755), and there is widespread but sparse low-level instrumental seismicity across the continent from the British Isles to Karelia, and paleoseismological works suggest several areas of active deformation (e.g., along the Rhine Graben (e.g., Camelbeeck et al., 2007; Grützner et al., 2016; van Balen et al., 2019), Lower Saxony Basin (e.g., Brandes et al., 2012; Brandes and Winsemann, 2013; Brandes et al., 2018; Muller et al., 2021), Cheb Basin (e.g., Štěpančíková et al., 2019), and the Sudetic Marginal Front (e.g., Štěpančíková et al., 2012; Štěpančíková et al., 2022)).

Similarly, North America, east of the Rocky Mountains and Cascades, is considered as a stable continental interior, largely seismically quiescent. However, there are a few notable areas of localized seismicity (e.g., the New Madrid Seismic Zone, the East Tennessee Seismic Zone, the St. Lawrence Valley Seismic Zone), although none of these have detectable ongoing tectonic strain accumulation associated with them (Craig and Calais, 2014; Kreemer et al., 2014; Boyd et al., 2015; Kreemer et al., 2018). The Teton and Yellowstone ranges in the central United States have previously been identified as hosting Holocene fault scarps with slip rates which correlate with modelled strain rates from variations in

local ice sheets (Hampel et al., 2007; Hampel et al., 2021).

In this work, we seek to quantify the time-dependent strain and stress rates in continental interiors associated with the evolution of the volume of the major northern hemisphere ice sheets, and how this may impact fault activation in Europe and North America. Our calculations focus on the European ice sheets (principally those over Fennoscandia, the Alps and the British Isles - see Fig. 1a) over ~40 ka, and the Laurentian ice-sheet of North America (see Fig. 3a).

Several studies have indeed suggested that the distal effects of the Fennoscandian deglaciation influenced fault behaviour of central Europe in the Holocene – Late Pleistocene. Houtgast et al. (2005) used variations in sedimentation rate across the Geelen Fault (Netherlands) to infer an increased slip-rate between 10 and 15 ka that they relate to glacially-induced variations in the regional deformation rate and related increase in fault activity, with nearby faults experiencing moderate-magnitude earthquakes during the same time period (van Balen et al., 2019). In northern Germany, the reactivation of faults in the Lower Saxony Basin, interpreted from the deformation of Pleistocene sediments, has been suggested to result from the development and decay of the Fennoscandian forebulge (Brandes et al., 2012; Brandes and Winsemann, 2013; Brandes et al., 2015; Mueller et al., 2020; Muller et al., 2021). In western Poland and the Czech Republic, recent work on the Sudetic Marginal Front, ~ 150 km from the maximum ice margin, indicates that this fault experienced increased slip rates during periods of glacial loading (Štěpančíková et al., 2022).

In North America, fewer studies have considered the interaction of ice sheets on fault systems, but examples do include New Madrid (Grollimund and Zoback, 2001), the Teton Ranges and Basin and Range (e.g., Hampel et al., 2007; Hampel et al., 2009; Hampel et al., 2010; Hampel et al., 2021) and Alaska (Sauber and Molnia, 2004; Sauber and Ruppert, 2008; Rollins et al., 2020; Sauber et al., 2021).

Here we will show that the far-field strain-rates resulting from changes in the ice load have been significantly greater in the past 25 ka than the slow rates of tectonic deformation currently taking place in continental Europe, and that they have migrated significantly over time. Whilst the mode of failure in earthquakes typically reflects the release of long-term tectonic stresses, and not the transient stresses induced by changing surface loads, their timing and location may be affected by these transients. Although the models presented here are non-unique, they provide quantitative estimates of strain and stress rate variations that should help in interpreting paleoseismic records for seismic hazard assessment where more detailed consideration of the role of non-tectonic processes has not yet been incorporated. This is particularly important for critical infrastructure – nuclear waste storage and disposal facilities, for instance – whose design is based on safety projections over very long time intervals (10³ to 10⁶ years), and which are typically sited in low-strain environments.

## 2. Modelling approach

To assess the effect of the redistribution of ice masses on continental strain rates in Europe and North America, we construct a series of models that allow us to calculate stress and strain that result from changes in surface loading over a glacial cycle, similar to the approach described in Craig et al. (2016) and Caron et al. (2017). Models are constructed under the assumption that the Earth behaves as a self-gravitating visco-elastic sphere (radius 6371 km). We calculate the response of the crust and mantle to a periodic surface load, expressed up to a spatial resolution of spherical harmonic degree 128, equating to a lateral resolution of ~300 km at the Earth's surface. Boundary conditions are specified at the core-mantle boundary (2891 km depth) and at the free surface, where changes in surface load are applied as a predetermined time-variable radial stress.

Unlike commonly used methods based on the computation of normal modes, our method is based on the Fourier decomposition of the time-
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Strain-rate distribution across Europe. (a) Ice volume at 20 ka from ANU-ICE. Solid contours are at 200 m intervals. Dashed contour is the 100 m contour, as a proxy for the ice margin. (b)-(f) Second invariant of the deviatoric strain-rate tensor at (b) 41–40 ka, (c) 31–30 ka, (d) 21–20 ka, (e) 11–10 ka, (f) 1–0 ka. This effectively shows the magnitude of the overall strain-rate. The scale used is the same in each case. All results are calculated at the free surface.</figcaption></figure>
dependent variation for each spherical harmonic component of the load. The response of the Earth for each spherical harmonic and each time-frequency is then computed using the classical method used for computing elastic Love numbers (Alterman et al., 1959; Cathles, 1975) except that the elastic parameters are replaced by complex numbers which represent the viscoelastic parameters as a function of frequency.

We use the ANU-ICE model (Lambeck et al., 2014) for changes in the extent and volume of major ice sheets through time. This ice model and our modelling approach are global in extent. We resampled the initial ice model onto a 1° × 1° spatial grid and to 1 ka time intervals, by linear interpolation. Since our modelling approach requires, for mathematical simplicity, that the surface load variation over the timescale of the model be periodic, a 250 ka loading cycle is supplemented by an additional 200 ka of no load change from the present, in order to allow for relaxation of the glacial process. Then the loading cycles are merged back into the time re-initialisation of glaciation at 250 ka to create a periodic signal.

Accumulation of the Fennoscandian ice sheet takes place over the late Pleistocene to the last glacial maximum at 23–20 ka. Then ice retreat takes place gradually until 10 ka, at which point deglaciation of

Fennoscandia is complete. In the British Isles, ice is concentrated over Scotland and areas of northern England, northern Ireland and Wales. It is connected to the main Fennoscandian ice sheet during peak glaciation, but with both the peak and final termination of major glaciation occurring slightly earlier, at ~25 ka and ~15 ka respectively. The Alpine ice sheet, whilst much more minor in amplitude and extent than the previous two, is important for strain patterns in central Europe. It peaked between ~24 and ~10 ka, with a relatively rapid decline accomplished by ~7 ka. In North America, the Laurentian ice sheet covered much of Canada and the northmost USA over the Pleistocene, peaking at ~20ka, before a more gradual, steady decline and retreat until end glaciation at around ~6 ka.

The ice loading model is adapted to account for the conjugate changes in oceanic loading. At the resolution of our model, fully solving the sea-level equation would produce only minor variations in the strain and stress fields. We instead implement broad-scale changes in oceanic loading by redistributing uniformly across the oceans the ice load removed without modifying coastlines, whilst conserving the total equivalent water load at all time steps. We do not recalculate coastlines at each time interval, and so exclude from our model the flooding of

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Strain-rate time series across Europe. (a) Second invariant of the deviatoric strain-rate tensor at 13–12 ka. (b)-(k) Profiles of the principal axes of the horizontal strain-rate tensor through time at the locations shown on (a). Points are coloured to indicate the orientation (in azimuth clockwise from north) of each axis. Note that the strain-rate scale is different on each profile. All results are calculated at the free surface.</figcaption></figure>
shallow continental shelf regions like Irish Sea, North Sea, English Channel, and northernmost Adriatic and the effect this would have on the near-field stress and strain fields. The exception to this is the loading of the Black Sea, which we model as being unconnected to the global oceanic system prior to 7 ka. At 7 ka, the opening of the Bosphorus Strait leads to the integration of the Black Sea back into the global oceanic system. This only has a secondary effect (compared to global sea-level changes) on the strain and stress fields of Anatolia around 7 ka.

The flooding of the Black Sea produces a notable kink in the strain-rate profile for Anatolia at 7 ka, as shown on Fig. 2, and has been suggested to play a major role in the stress state of Anatolia, particularly around the North Anatolian Fault (Luttrell et al., 2007). However, given the relatively small contribution of the Black Sea to the total oceanic volume, this has minimal effects on more distal regions, with no discernible associated kink in strain rate present in profiles on Fig. 2 at greater distances from the Black Sea. Hence, whilst the precise timing and rate of this Black Sea flooding remains a topic of some debate (Ryan et al., 2003), variations of a few kips do not significantly alter our model results. For simplicity, shallow endorheic oceans such as the Caspian Sea, Lake Chad, etc. are assumed to be disconnected from the global ice/ocean system, and their load-evolution is not incorporated into our model.

Elastic properties are taken from the seismologically-derived one-dimensional Preliminary Reference Earth Model (Dziewonski and Anderson, 1981) for a spherically-symmetric Earth. The 1-dimensional viscosity ($\eta$) structure used is based on that of Zhao et al. (2012), which comes twinned with the ANU-ICE model which we are also using. It incorporates a 101 km-thick elastic lithosphere over an upper mantle with $\eta = 4.2 \times 10^{20}$ Pa s, a lower mantle with $\eta = 1.0 \times 10^{22}$ Pa s, and a transition between the two at 660 km below the free surface. Comparisons to models constructed using the same approach from the ICE-5G ice history model (Peltier, 2004) and the twinned VM5a viscosity structure (Peltier and Drummond, 2008) demonstrate that, whilst the finer details of the strain and stress field generated do differ, the large-scale features which are the concern of this paper are found in both Earth/ice model pairs (e.g., Steffen et al., 2019). These small-scale differences are smaller than other unquantified effects such as that of failing to incorporate the 3-dimensional structure of both the elastic lithosphere and the visco-elastic underlying mantle.

The most problematic issue in such calculations results from the relatively poorly constrained viscosity of the lower mantle. Observational constraints on the viscosity of the lower mantle are largely derived from long-wavelength glacial isostatic adjustment (GIA hereafter), and viscosity is determined in conjunction with long-wavelength ice load history (e.g., Peltier, 2004; Zhao et al., 2012). For the Laurentian icesheet in North America, this poses a particular problem, due to the sheer scale of the ice sheet at its maximum extent, and the paucity of geological and geomorphological data from the continental interior to constrain this. Here, where we are mainly concerned with the far-field effects of ice-loading beyond the edges of the ice margin, the longer-wavelength impact of lower mantle viscosity is a particular problem. To test the impact of uncertainties in lower-mantle viscosity on the induced intraplate strain fields we show for North America, we also run tests, assessing how much these strain fields vary if we change the lower-mantle viscosity, increasing or decreasing it by factors of 5 and 10 (see Section 4.3).

Model time increments are set to 1000 yrs, with the full strain and stress tensors computed at each time interval. Strain- and stress-rate tensors are calculated by differencing the solutions for displacement at adjacent time-steps prior to the calculation of strain and stress tensors. The results shown in Figs. 1, 2, and 4 are for the strains at the free surface, showing the 2nd invariant of the deviatoric strain-rate tensor (effectively the magnitude of shear strain), and hence are comparable to those measurable at the surface by geodesy or paleoseismology.

### 3. Time/space-variable strain-rates at continental scale

Our model results (Figs. 1 and 2) show that whilst present-day glaciation-induced strain rates in Europe are low outside of Fennoscandia ($< 5 \times 10^{-9}$ yr$^{-1}$), they were significantly greater over much of the Holocene and late Pleistocene than they are at present. In addition, model results show that the strain-rate field was spatially complex (Fig. 1) from 40 to about 10 ka, a result of the interplay between the slightly asynchronous evolution of the Fennoscandia/Russian Arctic, British Isles, and Alpine ice sheets (Fig. 1a) and the influence of oceanic volume changes. Similarly, horizontal strain rates in North America associated with the growth and decay of Laurentian ice sheet reach $\sim 10^{-7}$ yr$^{-1}$ near the ice margins themselves, and exceed $\sim 10^{-8}$ yr$^{-1}$ in the continental interior, extending to the Central United States – far in excess of anything observable at the present day at such latitudes (Calais et al., 2006; Kreemer et al., 2014; Kreemer et al., 2018).

Changes in surface load result in an immediate elastic response, which dominates the deformation field at short-wavelengths, followed by a slower long-wavelength viscous response, the amplitude of which decays over time as the system re-equilibrates. Ongoing long-wavelength deformation at present in Fennoscandia and northern North America, some 10 ka after the end of major glaciation, is driven by this viscous response (Figs. 1f, 3c). The shorter-wavelength ice load over the Alps, for example, is instead predominantly supported elastically, and so produces a rapid, more localised solid-Earth response (Fig. 1e), with a smaller, delayed, viscous component.

Whilst the large-scale pattern of deformation shown on Fig. 1 may appear, to first-order, similar through time, Fig. 2 shows that the magnitude and orientation of the principal axes of the horizontal strain-rate tensor go through a number of rotations and reversals throughout the glacial cycle around the periphery of the major ice sheets. These reversals are most simply observed by considering central Turkey (Fig. 2k), a location far enough away from the major ice sheets that the model strain-rates are dominated by the effect of changing sea level in the Black Sea and the eastern Mediterranean rather than by variations of the continental ice mass. One of the principal axes of the horizontal strain-rate tensor is hence always oriented approximately east-west, with a low magnitude. The other axis is consistently oriented approximately north-south, but reverses from compression (positive values on Fig. 1)) to extension (negative values on Fig. 1) at around 19 ka, when the global continental ice mass transitions from increase to decrease, with a concomitant shift from sea-level fall to sea-level rise. The notable kink in the N/S-orientated axis at $\sim 7$ ka is due to the connection of the Black sea to the global ocean system, as previously discussed.

Peak strain-rates at any time-step correspond to the location of the largest changes in the surface load as they result from the immediate elastic and initial rapid viscous Earth response. Hence, the largest signal in Fig. 1c,d,e is observed within Fennoscandia, at the location of contemporaneous ice load change, and on Fig. 3b,c in the areas of Arctic Canada associated with the greatest thickness of the Laurentian icesheet. However, significant strain-rates reach far beyond the ice margins, with a long wavelength viscous response driving crustal deformation across central Europe and western Russia, and extending as far as the Balkans and the north Caspian basin. This large-scale viscous response persists long after the eventual decay of the ice load (Fig. 2).

Outside the ice margin, the most rapid strain-rate changes are produced instead by the growth and then decay of the Fennoscandian ice-sheet forebulge, where deformation is dominated by the elastic support of the ice margin lithosphere. This is best shown on Fig. 1b by the annular structure around the Norwegian coast, through the Baltic states and down to northern Poland, and on Fig. 1d by the sharp spike in strain rates through Eastern Europe and Karelia. For North America, this is most apparent on Fig. 3c, where the band of high-rate deformation that broadly aligns with the Canada/United States border reflects the ongoing collapse of the Laurentian forebulge – a feature detectable with
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Strain rates in North America. (a) Ice extents at 20 ka from ANU-ICE. Solid contours are at 500 m intervals. Lower panels show profiles as indicated on the map, with ice sheet thicknesses at 2 ka intervals from 20 ka to 6 ka. (b) Second invariant of the strain-rate tensor at 19–18 ka. Lower panels show profiles, indicating the magnitudes and (as symbol colour) orientation of the principal axes of the strain rate tensor. (c) as in (b), but for 1–0 ka. Note that the scale in the lower panels is reduced by 1/5 in comparison to (b).</figcaption></figure>
modern GNSS geodesy (e.g., Calais et al., 2006; Kreemer et al., 2018).

The growth and decay of this forebulge and the migration of the strain rate peak with ice growth and removal are particularly relevant to the time-variable strain-rates of both continental Europe and intraplate North America. In Russian Karelia (Fig. 2j), a brief period of rapid NW-SE extension between 24 and 19 ka, coincident with the development of the closest part of the Fennoscandian ice sheet at the LGM, is followed by a long interval of low-rate compression, reflecting the gradual decline of

ice along the northeastern margin of the ice sheet. A similar time-evolution is seen for the North Sea (Fig. 2b). In both of these locations within the Fennoscandian forebulge, model strain-rates are in excess of $5 \times 10^{-8} \text{ yr}^{-1}$, a value that would be easily measured using today's space geodetic techniques.

Across the rest of continental Europe, model strain-rates show significant variations in magnitude and orientation through time that may not be intuitive. In the northern Czech Republic, for example, in addition

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Strain-rate evolution in the European Cenozoic Rift System. (a) The Cenozoic European Rift System. Grey dots are earthquakes form the European-Mediterranean Earthquake Catalogue for 1000–2006 (Grünthal and Wahlström, 2012), filtered for $M_W > 3.5$, and scaled by magnitude. Black lines are the fault systems of the Upper and Lower Rhine Graben after (Vanneste et al., 2013), and the North German Basin after (Brandes et al., 2012). The sense of motion shown is based on the Cenozoic motion of the fault, and may differ from the sense of motion in recent earthquakes, where reactivation has occurred. Bel: Belgium. Nld: Netherlands. LRG: Lower Rhine Graben. URG: Upper Rhine Graben. LSB: Lower Saxony Basin. (b) - (f) Principal axes of the horizontal strain-rate tensor (coloured bars, blue for extension, red for compression), overlain on the second invariant of the deviatoric strain-rate tensor. The time interval displayed is shown in the top left corner of each panel. The scale for strain-rate crosses is multiplied by a factor of 2 on panel (e) and a factor of 5 on panel (f), to make the results visible. (g), (h), (i) Evolution of the principal axes of the horizontal strain-rate tensor for the Lower Rhine Graben, Upper Rhine Graben and Lower Saxony Basin, respectively. Point colour on (g), (h), (i) indicates the angle between the principal strain axis and each fault system.</figcaption></figure>
to variability in the strain-rate magnitude, model results also shows 45° rotation in the orientation of the tensor in $\leq 6$ kyrs (Fig. 2i). Similarly, Germany, within the forebulge of the Fennoscandian ice sheet and close enough to the Alps to be affected by the effects of Alpine glaciation, presents a complex evolution through time – discussed in more detail in Sections 4.1 & 4.2.

The effect of ocean margin loading is particularly visible along the coast of North Africa (Figs. 1c and 1e). This feature is dominated by the short-wavelength flexure of the margin, resulting in margin-perpendicular extension onshore and compression offshore during times of increasing oceanic volume (continental ice loss – e.g., Fig. 1e), and the converse during times of ocean volume decrease (continental ice accrual – e.g., Fig. 1c). The flexural effects of ocean margin loading, particularly with respect to strike-slip fault systems, has been previously investigated in detail elsewhere (e.g. Luttrell and Sandwell, 2010; Brothers et al., 2013). Whilst our modelling approach has a more limited spatial resolution and a more simplistic implementation of coastal loading in comparison with that of Luttrell and Sandwell (2010), ours has the advantage that we include long-wavelength effects due to the large-scale ice loads – necessary for regions within ~2000 km of the ice margin. In summary, Figs. 1–3 show that strain-rates induced by variations of continental ice masses are heterogeneous in both space and time in regions outside the ice margin. In addition, model results show that this process can result in strain-rates in these regions that are significantly larger than typical tectonic values in stable continental regions ($< 1 \times 10^{-9} \text{ yr}^{-1}$, Nocquet (2012), Calais et al. (2016)), reaching up to $20 \times 10^{-9} \text{ yr}^{-1}$ at the 1000-yr resolution of our model.

## 4. Regional examples

Although the above description of model results focuses on strain-rates, the activation of faults should more properly be discussed in terms of the stress, or the changes in stress, acting on them. However, correctly doing so requires a priori knowledge of the geometry and slip direction of faults in a given region, information that is rarely available in low-strain rate environments. Additionally, a robust test of the extent to which ice sheet load variations may modulate seismicity would require confronting modelling results with a complete paleoseismic catalogue spanning a period longer than the glacial cycle. Again, such an exhaustive paleoseismic catalogue is not yet available for either Europe or North America as a whole. In the following, we therefore focus on three of the best-studied areas of intraplate seismicity within continental Europe and North America in terms of paleoseismicity, the European Cenozoic Rift System (ECRS), the Lower Saxony Basin (LSB; Fig. 4a), and the New Madrid Seismic Zone (NMSZ; Fig. 6a). In all cases, significant effort has been put into establishing a paleoseismic record over the Holocene as well as the geometry and slip direction of the major potentially seismogenic faults (e.g., Kockel, 2003; Vanneste et al., 2013; Tuttle et al., 2005). We note that there are other regions within central and Northern Europe suggested to have been active over the Holocene (e.g., the Sorgenfrei-Tornquist zone, Brandes et al. (2015), Brandes et al. (2018)), but we focus on the ECRS, LSB, and NMSZ, where the fault dip and kinematics are both well known, and consistent across the fault system.

### 4.1. The European Cenozoic Rift System

The ECRS system stretches from the northern edge of the Alpine orogeny to the North Sea (Fig. 4a). It is split into two sections, the NNE-SSW trending Upper Rhine Graben (URG) and the NW-SE trending Lower Rhine Graben (LRG, also known as the Roer or Rur Valley Graben). The ECRS is one of the most seismically active areas of intraplate Europe and has been the locus of damaging earthquakes, including the $M_{L1}6.4$, 1756, Düren earthquake, the $M_{L1}5.8$, 1951, Euskirchen earthquake, and more recently the $M_{L1}5.1$, 1992, Roermond earthquake (Hinzen and Oemisch, 2001) with a damage cost estimated at 125

million euros. Seismic hazard within the ECRS is therefore of concern to a number of European nations, given the proximity of several major urban centres, including Strasbourg, Düsseldorf, Köln, and Eindhoven.

Geodetic measurements have so far not been able to detect significant tectonic strain across the ECRS (e.g. Nocquet, 2012; Fuhrmann et al., 2015), consistent with the low paleoseismic estimates of average Quaternary fault slip rates ($\leq 0.1 \text{ mm yr}^{-1}$, Vanneste et al. (2013)). Geologically-derived estimates for large earthquake recurrence intervals in the LRG range from 6 ka to $\geq 80$ ka (Vanneste et al., 2001; Vanneste et al., 2013; Grützner et al., 2016), and hence are comparable to, or longer than, the typical duration of a given orientation of the strain-rates shown in Figs. 1 and 2. Paleoseismic studies in the URG are more sparse, but indicate similar rates of motion of ~0.1–0.2 mm/yr (e.g., Meghaoui et al., 2001; Becker et al., 2005).

The LRG lies within the forebulge area of the Fennoscandian ice sheet (Fig. 2a,g), where model results show a transient episode of conglacial extension and deglaciation compression as the ice advances and retreats. The URG is also affected by the time-varying Fennoscandian ice load, but is close enough to the shorter-wavelength Alpine ice load that this has a additional effect. In addition, strain rates in the URG are likely affected by the ongoing erosion taking place across the Alpine orogenic belt, which produces a measurable geodetic strain signal (Sternai et al., 2019), but is not incorporated in our model.

Fig. 4 shows a close-up of the evolution of strain-rate in north-central Europe as a result of GIA over the past 25 ka. In order to determine whether GIA promotes fault activation of the ECRS bounding faults, we assume, to first order, that failure is promoted when one of principal strain-rate axes is both perpendicular to the fault orientation (points shaded black on the lower panels of Fig. 4) and is significantly negative, indicating an increase of the extensional strain.

We observe, for both the LRG and URG, a rather complex evolution of the principal axes of the strain rate tensor. At no point do our models indicate that these structures are subjected to simple rift-perpendicular extension. The three-dimensional nature of the strain-rate field rarely produces a strain-rate tensor consistent with uni-directional extension or compression. Even at times where one of the principal axes of the horizontal strain-rate tensor is negative and rift perpendicular, the other axis is typically positive to a similar magnitude and rift-parallel, as demonstrated for the LRG at 19–18 ka (Fig. 4c) and the URG over the last 1 ka (Fig. 4f).

In Fig. 5, we calculate rates of change in normal, shear, and Coulomb Failure stress on the LRG, URG, and LSB. All rifts are assumed to comprise pure-dip-slip normal faulting, at a dip of 60°. Coulomb Failure stresses are calculated using an effective coefficient of friction of 0.4. In terms of GIA-induced stress on rift-bounding faults, Fig. 5 indicates significantly larger temporal variations in the LRG than in the URG, predominantly due to its closer proximity to the Fennoscandian icesheet. Both grabens show time intervals where failure is enhanced or inhibited by the effects of GIA. In the URG, positive Coulomb stress changes never exceed 0.1 kPA/yr, indicating that the process modelled here likely had minimal impact on fault activation. In the LRG, increased hangingwall sedimentation rates from 15–10 ka have been suggested to be a result of an increase in fault activity (slip rate) during this time period due to the time-variable influence of post-glacial processes (Houtgast et al., 2005; van Balen et al., 2019). However, model Coulomb stress changes during this time interval show a (slight) decrease that does not support an increase in normal-faulting activity. Time intervals of increased model Coulomb stress, e.g., from 20–14 ka in the case of the LRG, are not correlated with documented enhanced fault activity, although we note that it corresponds to the reported age of the most recent earthquake on the Geelen fault in the LRG (Vanderberghe et al., 2009).

### 4.2. Lower Saxony Basin

The Lower Saxony Basin in northern Germany (LSB; Fig. 4), bounded by WNW-ESE trending faults, initially formed during the Permian as an
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Stressing-rate evolution in the European Cenozoic Rift System. Each panel shows the time-variation in glacially-induced stressing rate in terms of normal, shear, and a Coulomb Failure stress, for the Lower Rhine Graben (top panel), Upper Rhine Graben (middle panel), and Lower Saxony Basin (bottom panel). Stress is calculated at 10 km depth assuming planar faults with a geometry based on their surface strike, a dip angle of 60°, and pure dip-slip, normal faulting, motion.</figcaption></figure>

extensional rift system. Many of these faults were then reactivated as compressional thrust faults during basin inversion in the late Cretaceous-Paleocene (Kockel, 2003), most prominently the Osning thrust at the southern margin of the basin.

Trenching across the Osning thrust suggests that a more rapid interval of small-scale extension and inversion occurred over the last glacial cycle (Brandes et al., 2012; Brandes and Winsemann, 2013; Brandes et al., 2018), with a small amount of extensional slip on the fault during ice advance as the forebulge developed in northern Germany, followed by reversal and thrust motion on the same fault during and following deglaciation as the forebulge collapsed. Fig. 4 shows that LSB faults were indeed favourably aligned to the glacially-induced strain-rate field to undergo extension during ice accrual prior to ~20 ka, and then reversed to compression from about 16 – 8 ka. Model Coulomb stress changes on Fig. 5 are positive, hence consistent with fault activation, during the 16 – 8 ka time interval. However, this does not hold prior to ~20 ka. These results are similar to those of Brandes et al. (2015), who suggest that the removal of the Fennoscandian icesheet promoted reverse-faulting failure of the Osning thrust between 16 and 10 ka.

### 4.3. The New Madrid Seismic Zone

The New Madrid Seismic Zone (NMSZ; location on Fig. 3a) is a region of active intraplate seismicity within the continental interior of North

America. Whilst present-day seismicity is typically <M4, the area experienced a sequence of large-magnitude (M>7) earthquakes in the winter of 1811–1812 (Johnston, 1886; Hough et al., 2000), with geological evidence for other major earthquakes during the later Holo-cene (Gold et al., 2019), with up to six episodes of regional liquifaction since ~3000 BČE. (Tuttle et al., 2005; Holbrook et al., 2006). Present-day strain rates in the NMSZ are undetectable ~1–3 × 10⁻⁹ yr⁻¹ (Craig and Calais, 2014; Boyd et al., 2015), leaving the causes of this concentration of intraplate seismicity uncertain. Here, we do not attempt to answer this question, but instead use New Madrid as an example region to investigate the impact of far-field ice-loading on intraplate strain. In Fig. 6, we show time-series for strain- and stress-rates at New Madrid driven by GIA, and three snapshots of the strain field.

Several previous studies suggested or investigated the impact of GIA-related deformation on the NMSZ (e.g., Grollimund and Zoback, 2001; Hough and Page, 2011). Unlike Grollimund and Zoback (2001), we do not include a specific rheologically-weak zone beneath the NMSZ. In Grollimund and Zoback (2001), this serves to focus GIA-induced strain into the region of the NMSZ, producing strain rates capable of producing repetitive seismicity. We instead continue with the radially-symmetric rheological model as described in Section 2, focusing on the longer-wavelength impacts of GIA across the continental interior.

The NMSZ consists of a NE-striking, right-lateral strike-slip fault, and
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** The New Madrid Seismic Zone. (a) Seismicity in the New Madrid Seismic Zone, from the CERI catalogue. Red circles are the approximate locations of the 1811–1812 earthquakes. (b) – (d) show the second invariant of the strain-rate tensor (as shading), overlain by principal axes of the strain-rate tensor (blue for extension, red for compression), at three time intervals. (e) shows a time series for the magnitude and (as point colouring) orientation of the principal axes of the strain-rate tensor (expressed as azimuth) at the location of the NMSZ. (f) is as in (e), but for the stressing-rate tensor.</figcaption></figure>

a SW-dipping, SE-striking reverse fault, both of which likely ruptured in the 1811–1812 earthquake sequence. Interestingly, our modelling suggests that the strain and stress fields induced by changes in ice-loading in this region, although far too small to have loaded the faults sufficiently in and of themselves, would have been consistent with promoting failure of the strike-slip system between 18–6 ka, and then promoting failure of the reverse fault system from 5–0 ka, in keeping with paleoseismic evidence suggesting persistent failure over the later Holocene (e.g., Tuttle et al., 2005; Holbrook et al., 2006; see Fig. 6e). This contrasts with the earlier findings of Wu and Johnston (2000), who predicted the promotion of failure in the NMSZ to have only started from only 200 yrs ago, with the difference likely resulting from the difference in ice model and viscosity structure used, particularly in the lower mantle. Whilst other processes (tectonic or otherwise) must have been involved in loading the faults of the NMSZ to the stage of failure, and are required to explain why earthquakes are concentrated around the NMSZ, and not elsewhere in the continental interior, the removal of the Laurentian ice sheet, under the assumptions made here, would have moved the NMSZ closer to failure.

As discussed in Section 2, the deeper viscosity of the mantle plays a dominant role in controlling the longest-wavelengths of induced deformation. However, these viscosities remain poorly constrained, leading to significant uncertainty in the magnitude and decay timescale of the far-field GIA signal – particularly the horizontal components of the strain tensor. To rigorously test the impact that uncertainties in the lower mantle viscosity have on the surface deformation field, varying

the viscosity structure should be coupled with a re-determination of the ice history, as the two are derived in combination. Such an endeavour is beyond the scope of our study. Instead, as a test for the impact that uncertainties in lower mantle viscosity may have, we modify the lower mantle viscosity in the structure determined in Zhao et al. (2012), as detailed in Fig.S1. As this figure demonstrates, variations in lower mantle viscosity have a major impact on the magnitude of the principal axes of the horizontal strain-rate tensor, with much faster decay in far-field strain-rates for a reduced viscosity. However, the times at which changes are seen in the orientation of far-field strain-rates is more closely related to changes in the growth/decay rate of the ice load, and is relatively insensitive to viscosity.

## 5. Continental margin loading

The effect of changing ocean volumes as a result of variations in continental ice masses on near-marginal faulting has been studied previously, with a particular emphasis on near-coastal transform fault systems (Luttrell and Sandwell, 2010), and marginal fault-related margin slope failure (Brothers et al., 2013). However, changing ocean volumes, and the strain-fields induced by the resulting flexure of the margin, may affect a wide range of active near-margin fault systems. As shown on Fig. 1, strain-rate variations induced by this process can be observed in the model results for the tectonically-active regions of the Atlas margin in North Africa, and the N-S orientated extensional system of western Anatolia.
For instance, Fig. 2.k illustrates the strain-rate evolution at the eastern end of the extensional systems of Anatolia, in central Turkey. There, the ocean-induced strain field is dominated by the flexure of Anatolia as the volumes of the Black Sea and Eastern Mediterranean vary. Model calculations show little variation in E-W strain, but N-S strain-rates that vary between $\pm 5 \times 10^{-9} \text{ yr}^{-1}$. As the geodetically observed present-day strain-rates in that same area are estimated to be around $25 \times 10^{-9} \text{ yr}^{-1}$ (Nocquet, 2012; Piña-Valdés et al., 2022), ocean loading-induced strain may lead to fluctuations of about 20% of the overall extension rates. As a result, one may expect increased rates of seismicity during times when oceanic loading leads to N-S extension, in agreement with the regional tectonics (e.g., 18–7 ka), and decreased earthquake occurrence when the opposite is the case (e.g., 29–20 ka).

Similar magnitudes of ocean-loading derived strain-rate are predicted for other active areas, such as Central Greece and peninsular Italy. However, their effect on seismicity rates is likely to be much smaller, due to the significantly greater tectonic strain-rates, in some cases exceeding $100 \times 10^{-9} \text{ yr}^{-1}$ (Nocquet, 2012; Piña-Valdés et al., 2022), and due to less favourable alignments between the secondary and tectonic strain fields than seen in western Anatolia.

An alternative example arises from considering the margins of North Africa through Morocco, Algeria, and Tunisia. In these regions – too distal from the major ice sheets for much of a direct deformation signal from changes in glacial loading – the major source of deformation is the elastic deformation associated with the changing water levels in the Mediterranean. As such, a simplistic load-induced stress field emerges (visible on Fig. 1c,e, in particular), in which, as water level rises, the onshore areas will be subject to an N-S extensional shallow stress change, with deeper N-S compression, which reverse during times of sea level fall. As these regions of North Africa are tectonically active, these induced stress fields, although likely small in comparison to the tectonic stresses, may have a minor modulating effect of the stress accumulation of faults in the region.

The values and wavelengths of the deformation associated with continental margin loading found here are however dependent on the shallow rheological structure, which is not accounted for in the global model used here. As the model parameters used here depend on fitting large-scale observations of glacial isostatic adjustment over continental ice masses that largely coincide with cratonic areas (e.g., Zhao et al., 2012), its average rheology is likely to be stronger, at lithospheric depths, than the non-cratonic continental margins described above. To fully understand the influence of both distal icesheet variations and ocean-loading requires more complex modelling, incorporating regional (and regionally-variable) rheological structures, and, particularly for the ocean-loading problem, the full solution of the sea level equation with time-variable coastlines and topography (Gomez et al., 2018; Whitehouse et al., 2019).

# 6. Implications for the ‘seismic cycle’

Seismic hazard assessment in continental interiors is often predicated on the assumption that faults behave in a quasi-steady-state manner in which they (1) accumulate stress over time at a steady rate dictated by long-term tectonics, then (2) release the accumulated stress in an earthquake when the shear stress on the fault exceeds its failure limit. In such a model, and in the absence of significant forcing other than long-term tectonics, seismic hazard can therefore be addressed by estimating fault slip rate from space geodesy or paleoseismology and extrapolating it to an earthquake recurrence time and/or an estimated earthquake population (e.g., Rollins and Avouac, 2019; Gerstenberger et al., 2020).

We have shown that strain – and hence for an elastic material, stressing – rates likely varied significantly in time and space in continental interiors as a result of glacial isostatic adjustment accompanying variations in icesheet volumes. For instance, in the three cases shown in

Fig. 5, significant GIA-related strain-rate variations between 40 and 10 ka are followed by negligible variations from ~10 ka onward. Hence, seismicity rates in the late Pleistocene and the Holocene may not necessarily be similar to each other for the same fault system. More generally, in areas where non-tectonic processes such as GIA cause significant time-variable strain-rates, the extrapolation of observational, historical, or paleoseismic data – the latter two usually being limited in terms of the number of earthquakes considered – to the present-day seismic hazard comes with the risk of mis-representing which faults are truly active tectonic structures, without additional consideration of what other non-tectonic processes may be impacting on regional earthquake occurrence. This may lead to the overestimate of hazard associated with faults with paleoearthquakes linked to transient processes like GIA, and underestimation of hazards now experiencing a more recent transient in regional deformation.

The magnitude of stress and strain rates induced by GIA are small compared to tectonic strain rates at plate boundaries or even in slowly deforming regions (typically well in excess of $10^{-8} \text{ yr}^{-1}$; Kreemer et al. (2014)). Moreover, the resulting strain and stress regime can alternate between compression, extension, or strike-slip over short time intervals (Fig. 2). It is therefore unlikely that GIA stresses by themselves can bring a fault to its point of failure, especially at distance from the principal load. However, if most crustal faults are in a state of failure equilibrium and if elastic strain is stored in the bulk of crust (e.g., Zoback and Healy, 1992; Townend and Zoback, 2000), including in stable continental interiors (Craig et al., 2016), then small stress perturbations caused by GIA may be sufficient to modulate and/or trigger seismicity. The stress changes involved are indeed similar to time-dependent stresses caused by hydrological loading that have been demonstrated to modulate seismicity in a variety of tectonic contexts, including stable continental interiors (e.g., Bollinger et al., 2007; Christiansen et al., 2007; Johnson et al., 2017; Craig et al., 2017; Rollins et al., 2020; Hsu et al., 2021).

Fig. 7 illustrates in a schematic manner how the superposition of a time-variable and a linear background tectonic stressing-rate may affect the timing of earthquake occurrence in a given area. We assume that earthquakes repeat for the same amount of accumulated stress within a given area and that there is always a favourably oriented fault able to rupture when that state is reached. The total stress build-up is the sum of the time-variable stressing-rate and of a linear, background, tectonic stressing-rate. The latter may be extremely small in stable continental regions, where strain rates are typically $< 2 \times 10^{-9} \text{ yr}^{-1}$ (Kreemer et al., 2018; Ding et al., 2019; Masson et al., 2019).

This simple conceptual model has several corollaries:

* Firstly, the presence of time-dependent stress obviously advances or delays the occurrence of earthquakes compared to a model where only tectonic stress is acting. This introduces a variability in the inter-event time compared to a theoretical, purely steady-state, system, in which earthquake occurrence would be regular and monotonic.

* Secondly, the variability of the inter-event time depends on the amplitude of the time-dependent stress changes with respect to the constant background tectonic stressing-rate. At the limit, if the latter is extremely small, such as in stable continental regions, then inter-event time depends solely on non-tectonic, time-dependent stress changes and may be very variable, and potentially non-repetitive. Conversely, if the tectonic stressing rate is large compared to time-dependent stress changes, such as at an active plate boundary, inter-event times will be much less variable as they are mostly dictated by the background tectonic loading. Whilst exaggerated to illustrate the point, in the simple example shown in Fig. 7, the inter-event time varies by ~50% – real values are doubtless substantially smaller.

* Thirdly, the superposition of the time-variable signal results in time intervals where the failure of well-oriented faults may be promoted (advanced) or delayed. In cases where the amplitude of the time-
<figure><img src="image_10.png" /></figure>

**Fig. 7.** Schematic stress accumulation in continental in-teriors. Simple model for the combination of a uniform background ‘tectonic’ stressing rate, and a superimposed time-variable ‘non-tectonic’ stressing rate. Green lines are for a time-invariant tectonic stressing rate, blue for a time-variable stressing rate, and red for the combined stress as seen by the fault. On the upper panel, turquoise lines indicate earthquakes (shown by black stars), assumed to occur at repeats of the same accumulated total stress, but which occur at variable intervals in model time. On the lower panel, grey-shaded regions indicate time periods where the combined stressing rate is negative, indicating that the fault is unlikely to rupture during these periods, despite the tectonic stress field.

variable signal exceeds that of the background stressing rate, this can go so far as to produce time intervals where the fault failure is inhibited.

Note that the illustrative model shown here in Fig. 7 treats failure as a simple threshold process, and includes no complex fault mechanics. The periods involved are long enough that processes relating to the nucleation of individual earthquakes are unlikely to matter. However, the frictional processes governing the accumulation, maintenance, and release of stress on individual fault planes are likely to lead to further complexity and variability in the temporal distribution of earthquakes on faults where such secondary processes are present that we do not attempt to quantify here. In the particular case of glacially-related load changes, there are also potential issues relating to fluxing of glacially-derived fluids through the upper crust, and the resulting changes in pore-fluid pressures, that we also do not consider in our simple model.

In the context of the GIA-related strain- and stress-fields that have formed the basis of this study, we illustrate that the role these processes may play on modulating seismicity and seismicity-rates may extend far

beyond the regions typically considered to be subject to GIA. Recent revision of the criteria for ‘Glacially-triggered faulting’ by Steffen et al. (2021) recognised that relevant areas may extend “several hundred kilometres” beyond formerly glaciated areas. However, as we demonstrate here, surface loading can impact crustal stresses and strains at distances up to ~1 wavelength from the load, with a resulting impact on the potential occurrence of seismicity. In the case of a globally-connected system (such as GIA) full consideration requires a truly global approach, incorporating all potential sources of stress (see Section 5).

## 7. Conclusion

We have demonstrated how strain-rates vary in space and time in Europe and North America solely as a result of the growth and decay of the Eurasian and Laurentian ice sheets since 40 ka. We show that such non-tectonic forcing can significantly influence the overall strain-rate field, and hence stresses that apply on faults within roughly one wavelength of the ice margin, in a rather complicated manner that includes
both the effects of changes in ice and ocean mass distributions.

Overall, the time-dependent pattern of GIA-induced strain-rate variations in Europe is dominated by the variability of the mass of the Fennoscandian icesheet, with smaller contributions from British Isles and Alpine glaciers. Continental margin loading as a result of icesheet melting adds a secondary complexity to the strain-rate variation pattern. Deformation comprises both the immediate elastic response to changes in load, particularly dominant at short wavelengths, and the viscous response, which dominates at longer wavelengths and over longer timescales. Model results indicate that strain-rates – and hence stresses that apply on faults – can be significant, with large spatial and temporal variations, during the late Pleistocene and peaking around the time of LGM. In some cases, the induced crustal stressing rates likely exceed the local tectonic stressing rates. Variations are much smaller over the Holocene, with the decay of major postglacial deformation across Europe, and are generally negligible after about 6 ka.

In regions where the background tectonic stressing rates are similar to, or smaller than, the superimposed non-tectonic rates, such effects can lead to time intervals where fault failure is advanced, delayed, or inhibited, depending on the alignment of the given fault system with the overall stress field. As a result, earthquake occurrence within given fault systems may become irregular, with long intervals of quiescence or bursts of enhanced activity. Whilst we lack sufficient paleoseismological data for a full assessment of the degree to which such variations influenced seismicity over this period, we recommend consideration of such effects in low-strain environments, as they add an additional uncertainty when using either modern-day geodetic strain rate fields, seismological records, or paleoseismic slip-rates based on small numbers of earthquakes, for long term seismic hazard assessment.

## CRediT authorship contribution statement

T.J. Craig: Conceptualization, Formal-analysis, Funding-acquisition, Investigation, Methodology, Validation, Visualization, Writing-original-draft. E. Calais: Conceptualization, Funding-acquisition, Supervision, Writing-review-editing. L. Fleitout: Conceptualization, Funding-acquisition, Methodology, Software, Writing-review-editing. L. Bollinger: Conceptualization, Funding-acquisition, Writing-review-editing. O. Scotti: Conceptualization, Funding-acquisition, Writing-review-editing.

## 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.

## Data availability

No data was used for the research described in the article.

## Acknowledgements

This work was funded through the French Investment Program SINAPS@ project by the Commisariat à l'Énergie Atomique and the Institute de Radioprotection et Sûreté Nucléaire, and was hosted by the LRC Yves Rocard (Laboratoire de Recherche Conventionné CEA-ENS-CNRS). TJC also thanks the Royal Society (under URF\R1\180088) for financial support during the final stages of this project. EC acknowledges funding from the Institut Universitaire de France. We thank K. Lambeck for making the ANU-ICE model available, and C. Gruetzner for his assistance with Rhine Graben data. Figures were created using the Generic Mapping Tools software package. We thank the editor and four anonymous reviewers for their comments, which have helped improve the manuscript.

# Appendix A. Supplementary data

Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.tecto.2023.229815.

## References

Alterman, Z., Jarosch, H., Pekker, C.L., 1959. Osciallations of the Earth. Proc. R. Soc. A 252. https://doi.org/10.1098/rspa.1959.0138.

Amos, C.B., Audet, P., Hammond, W.C., Bürgmann, R., Johanson, I.A., Blewitt, G., 2014. Uplift and seismicity driven by groundwater depletion in central California. Nature 509, 483–486. https://doi.org/10.1038/nature13275.

Becker, A., Ferry, M., Monecke, K., Schnellmann, M., Giardini, D., 2005. Multiarchive paleoseismic record of late Pleistocene and Holocene strong earthquakes in Switzerland. Tectonophysics 400. https://doi.org/10.1016/j.tecto.2005.03.001.

Bettinelli, P., Avouac, J.-P., Flouzat, M., Bollinger, L., Ramillien, G., Rajaure, S., Sapkota, S., 2007. Seasonal variations of seismicity and geodetic strain in the Himalaya induced by surface hydrology. Earth Planet. Sci. Lett. 266, 332–344. https://doi.org/10.1016/j.epsl.2007.11.021.

Bollinger, L., Perrier, F., Avouac, J.-P., Sapkota, S., Gautam, U., Tiwari, D.R., 2007. Seasonal modulation of seismicity in the Himalaya of Nepal. Geophys. Res. Lett. 34 https://doi.org/10.1029/2006GL029192.

Boyd, O.S., Smalley Jr., R., Zeng, Y., 2015. Crustal deformation in the New Madrid seismic zone and the role of postseismic processes. J. Geophys. Res. 120, 5782–5803. https://doi.org/10.1002/2015JB012049.

Brandes, C., Winsemann, J., 2013. Soft-sediment deformation structures in NW Germany caused by Late Pleistocene seismicity. Int. J. Earth Sci. 102, 2255–2274. https://doi.org/10.1007/s00531-013-0914-4.

Brandes, C., Winsemann, J., Roskosch, J., Meinsen, J., Tanner, D.C., Frechen, M., Steffen, H., Wu, P., 2012. Activity along the Osning Thrust in Central Europe during the Lateglacial: ice-sheet and lithosphere interactions. Quatern. Sci. Rev. 38, 49–62. https://doi.org/10.1016/j.quascirev.2012.01.021.

Brandes, C., Steffen, H., Steffen, R., Wu, P., 2015. Intraplate seismicity in northern Central Europe is induced by the last glaciation. Geology 43, 611–614. https://doi.org/10.1130/G36710.1.

Brandes, C., Steffen, H., Sandersen, P.B.E., Wu, P., Winsemann, J., 2018. Glacially induced faulting along the NW segment of the Sorgenfrei-Tornquist Zone, northern Denmark: Implication for neotectonics and Lateglacial fault-bound basin formation. Quatern. Sci. Rev. 189, 149–168. https://doi.org/10.1016/j.quascirev.2018.03.036.

Brothers, D.S., Luttrell, K.M., Chaytor, J.D., 2013. Sea-level-induced seismicity and submarine landslide occurrence. Geology 41 (9), 979–982.

Calais, E., Stein, S., 2009. Time-Variable Deformation in the New Madrid Seismic Zone. Science 323, 1442. https://doi.org/10.1226/science.1168122.

Calais, E., Mattioli, G., DeMets, C., Noquet, J.-M., Stein, S., Newman, A., Rydelek, P., 2005. Tectonics strain in plate interiors? Nature 438, E9–E10. https://doi.org/10.1038/nature04428.

Calais, E., Han, J.Y., DeMets, C., Noquet, J.M., 2006. Deformation of the North American plate interior from a decade of continuous GPS measurements. J. Geophys. Res. 111 https://doi.org/10.1029/2005JB004253.

Calais, E., Freed, A.M., Van Arsdale, R., Stein, S., 2010. Triggering of New Madrid Seismicity by late-Pleistocene erosion. Nature 466. https://doi.org/10.1038/nature09258.

Calais, E., Camelbeeck, T., Stein, S., Liu, M., Craig, T.J., 2016. A new paradigm for large earthquakes in stable continental plate interiors. Geophys. Res. Lett. 43, 10621–10637. https://doi.org/10.1002/2016GL070815.

Camelbeeck, T., Vanneste, K., Alexandre, P., Verbeeck, K., Petermans, T., Rosset, P., Evaraerts, M., Warnant, R., Van Camp, M., 2007. Relevance of active faulting and seismicity studies to assessment of long-term earthquake activity and maximum magnitude in intraplate northwest Europe, between the Lower Rhine Embayment and the North Sea. Geol. Soc. Am. Spec. Pap. 425, 193–224. https://doi.org/10.1130/2007.2425(14).

Caron, L., Métivier, L., Greff-Leffitz, M., Fleitout, L., Rouby, H., 2017. Inverting Glacial Isostatic Adjustment signal using Bayesian framework and two linearly relaxing rheologies. Geophys. J. Int. 209, 1126–1147. https://doi.org/10.1093/gji/ggx083.

Cathles, L., 1975. Viscosity of the Earths Mantle. Princeton University Press.

Christiansen, L.B., Hurwitz, S., Ingebritsen, S.E., 2007. Annual modulation of seismicity along the San Andreas Fault near Parkfield, CA. Geophys. Res. Lett. 34 https://doi.org/10.1029/2006GL028634.

Craig, T.J., Calais, E., 2014. Strain accumulation in the New Madrid and Wabash Valley Seismic Zones from 14 years of continuous GPS observation. J. Geophys. Res. 119, 1–20. https://doi.org/10.1002/2014JB011498.

Craig, T.J., Calais, E., Fleitout, L., Bollinger, L., Scotti, O., 2016. Evidence for the release of long-term tectonic strain stored in continental interiors through intraplate earthquakes. Geophys. Res. Lett. 43, 6826–6836. https://doi.org/10.1002/2016GL069359.

Craig, T.J., Chanard, K., Calais, E., 2017. Hydrologically-driven crustal stresses and seismicity in the New Madrid Seismic Zone. Nat. Commun. 8 https://doi.org/10.1038/s41467-017-01696-w.

Ding, K., Freymueller, J.T., He, P., Wang, Q., Xu, C., 2019. Glacial Isostatic Adjustment, Intraplate Strain and Relative Sea Level Changes in the Eastern United States. J. Geophys. Res. 124, 6056–6071. https://doi.org/10.1029/2018JB017060.

Dziewonski, A.M., Anderson, D.L., 1981. Preliminary reference Earth model. Phys. Earth Planet. Inter. 25 https://doi.org/10.1016/0031-9201(81)90046-7.
Fuhrmann, T., Caro Cuenca, M., Knöpfler, A., van Leijen, F.J., Mayer, M., Westerhaus, M., Hanssen, R.F., Heck, B., 2015. Estimation of small surface displacements in the Upper Rhine Graben area from a combined analysis of PS-InSAR, levelling and GNSS data. Geophys. J. Int. 203, 614–631. https://doi.org/10.1093/gji/ggv328.

Gerstenberger, M.C., Marzocchi, W., Allen, T., Pagani, M., Adams, J., Danciu, L., Field, E. H., Fujiwara, H., Luco, N., Ma, K.-F., Meletti, C., Petersen, M.D., 2020. Probabilistic Seismic Hazard Analysis at Regional and National Scales: State of the Art and Future Challenges. Rev. Geophys. 58 https://doi.org/10.1029/2019RG000653.

Gold, R.D., DuRoss, C.B., Delano, J.E., Jibson, R.W., Briggs, R.W., Mahan, S.A., Williams, R.A., Corbett, D.R., 2019. Four Major Holocene Earthquakes on the Reelfoot Fault Recorded by Sackungen in the New Madrid Seismic Zone, USA. J. Geophys. Res. 124, 3105–3126. https://doi.org/10.1029/2018JB016806.

Gomez, N., Latychev, K., Pollard, D., 2018. A Coupled Ice Sheet-Sea Level Model Incorporating 3D Earth Structure: Variations in Antarctica during the Last Deglacial Retreat. J. Clim. 31, 4041–4054. https://doi.org/10.1175/JCLI-D-17-0352.1.

Grollimund, B., Zoback, M.D., 2001. Did deglaciation trigger intraplate seismicity in the New Madrid seismic zone. Geology 29, 175–178.

Grünthal, G., Wahlström, R., 2012. The European-Mediterranean Earthquake Catalogue (EMEC) for the last millennium. J. Seismol. 16, 535–570. https://doi.org/10.1007/s10950-012-9302-y.

Grützner, C., Fischer, P., Reichertner, K., 2016. Holocene surface ruptures of the Rurrand Fault, Germany – insights from palaeoseismology, remote sensing and shallow geophysics. Geophys. J. Int. 204, 1662–1677. https://doi.org/10.1093/gji/ggv558.

Hampel, A., Hetzel, R., Densmore, A.L., 2007. Postglacial slip-rate increase on the Teton normal fault, northern region and Range Province, caused by melting of the Yellowstone ice cap and deglaciation of the Teton Range? Geology 35, 1107–1110. https://doi.org/10.1130/G24093A.1.

Hampel, A., Hetzel, R., Maniatis, G., Karow, T., 2009. Three-dimensional numerical modeling of slip rate variations on normal and thrust fault arrays during ice cap growth and melting. J. Geophys. Res. 114 https://doi.org/10.1029/2008JB006113.

Hampel, A., Karow, T., Maniatis, G., Hetzel, R., 2010. Slip rate variation on faults during glacial laoding and post-glacial unloading: implications for the viscosity structure of the lithosphere. J. Geol. Soc. London 167, 385–399. https://doi.org/10.1144/0016-76492008-137.

Hampel, A., Hetzel, R., Erdmann, M.-S., 2021. Postglacial slip distribution along the Teton normal fault (Wyoming, USA) derived from tectonically offset geomorphological features. Geosphere 17. https://doi.org/10.1130/GES02370.1.

Heki, K., 2003. Snow load and seasonal variation of earthquake occurrence in Japan. Earth Planet. Sci. Lett. 207, 159–164.

Hinzen, K.-G., Oemisch, M., 2001. Location and Magnitude from Seismic Intensity Data of Recent and Historic Earthquakes in the Northern Rhine Area, Central Europe. Bull. Seismol. Soc. Am. 91, 40 +56. https://doi.org/10.1785/0120000036.

Holbrook, J., Autin, W.J., Rittenour, T.M., Marshak, S., Goble, R.J., 2006. Stratigraphic evidence for millennial-scale temporal clustering of earthquakes on a continental-interior fault: holocene Mississippi River floodplain deposits, New Madrid seismic zone. USA. Tectonophysics 420 (3–4), 431–454.

Hough, S.E., Page, M., 2011. Toward a consistent model for strain accrual and release for the New Madrid Seismic Zone, central United States. J. Geophys. Res. 116 https://doi.org/10.1029/2010JB007783.

Hough, S.E., Armbruster, J.G., Seeber, L., Hough, J.F., 2000. On the modified Mercalli intensities and magnitudes of the 1811–1812 New Madrid earthquakes. J. Geophys. Res.: Solid Earth (1978–2012) 105 (B10), 23839–23864.

Houtgast, R.F., Van Balen, R.T., Kasse, C., 2005. Late Quaternary evolution of the Feldbiss Fault (Roer Valley Rift System, the Netherlands) based on trenching, and its potential relation to glacial unloading. Quatern. Sci. Rev. 24, 491–510. https://doi.org/10.1016/j.quascirev.2004.01.021.

Hsu, Y.-J., Kao, H., Burgmann, R., Lee, Y.-T., Huang, H.-H., Hsu, Y.-F., Wu, Y.-M., Zhuang, J., 2021. Synchronized and asynchronous modulation of seismicity by hydrological loading: A case study in Taiwai. Sci. Adv. 7 https://doi.org/10.1126/sciadv.abf7282.

Johnson, Christopher W., Yuning, Fu, Bürgmann, Roland, 2017. Seasonal water storage, stress modulation, and California seismicity. Science 336, 1161–1164.

Johnston, A.C., 1996. Seismic moment assessment of stable continental earthquakes – III. 1811–1812 New Madrid, 1886 Charleston, and 1755 Lisbon. Geophys. J. Int. 126, 314–344.

Karow, T., Hampel, A., 2010. Slip rate variations on faults in the Basin-and-Range Province caused by regression of Late Pleistocene Lake Bonneville and Lake Lahontan. Int. J. Earth Sci. 99, 1941–1953. https://doi.org/10.1007/s00531-009-0496-3.

Kierulf, H.P., Steffen, H., Simpson, M.J.R., Lidberg, M., Wu, P., Wang, H., 2014. A GPS velocity field for Fennoscandia and a consistent comparison to glacial isostatic adjustment models. J. Geophys. Res. 119 https://doi.org/10.1002/2013JB010889.

Kockel, F., 2003. Inversion structures in Central Europe – expressions and reasons, and open discussion. Neth. J. Geosci. 235, 277–291.

Kreemer, C., Blewitt, G., Klein, E.C., 2014. A geodetic plate motion and Global Strain Rate Model. Geochem. Geophys. Geosyst. 14, 3849–3889. https://doi.org/10.1002/2014GC005407.

Kreemer, C., Hammond, W.C., Blewitt, G., 2018. A Robust Estimation of the 3-D Intraplate Deformation fo the North American Plate From GPS. J. Geophys. Res. 123, 4388–4412. https://doi.org/10.1029/2017JB015257.

Lagerbäck, R., Sundh, M., 2008. Early Holocene faulting and paleoseismicity in Northern Sweden. Sveriges geologiska undersökning – Research paper, Uppsala, Sweden, C 836.

Lambeck, K., Rubby, H., Purcell, A., Sun, Y., Sambridge, M., 2014. Sea level and global ice volumes from the Last Glacial Maximum to the Holocene. Proc. Nat. Acad. Sci. 111, 15296–15303. https://doi.org/10.1073/pnas.1411762111.

Liu, Mian, Stein, Seth, 2016. Mid-continental earthquakes: Spatiotemporal occurrences, causes, and hazards. Earth Sci. Rev. 162, 364–386.

Luttrell, K., Sandwell, D., 2010. Ocean loading effects on stress at near shore plate boundary fault systems. J. Geophys. Res. 115 https://doi.org/10.1029/2009JB006541.

Luttrell, K., Sandwell, D., Smith-Konter, B., Bills, B., Bock, Y., 2007. Modulation of the earthquake cycle at the southern San Andreas fault by lake loading. J. Geophys. Res. 112 https://doi.org/10.1029/2006JB004752.

Masson, C., Mazzotti, S., Vernant, P., Doerflinger, E., 2019. Extracting small deformation beyond individual station precision from dense Global Navigation Satellite System (GNSS) networks in France and western Europe. Solid Earth 10, 1905–1920. https://doi.org/10.5194/de-10-1905-2019.

Mazzotti, S., James, T.S., Henton, J., Adams, J., 2005. GPs crustal strain, postglacial rebound, and seismic hazard in eastern North America: The Saint Lawrence valley example. J. Geophys. Res. 110 https://doi.org/10.1029/2004JB003590.

Meghraoui, M., Delouis, B., Ferry, M., Giardini, D., Huggenberger, P., Spottke, I., Granet, M., 2001. Active Normal Faulting in the Upper Rhine Graben and Paleoseismic Identification of the 1356 Basel Earthquake. Science 293, 2070–2073. https://doi.org/10.1126/science.1016018.

Müller, K., Polom, U., Winsemann, J., Steffen, H., Tsukamoto, S., Guenther, T., Igel, J., Spies, T., Lege, T., Frechen, M., Franzke, H.-J., Brandes, C., 2020. Structural style and neotectonic activity along the Harz Boundary Fault, northern Germany: a multimethod approach integrating geophysics, outcrop data and numerical simulations. Int. J. Earth Sci. 109, 1811–1835. https://doi.org/10.1007/s00534-020-01874-0.

Muir-Wood, R., 1989. Extraordinary deglaciation reverse faulting in northern Fennoscandia. In: Gregersen, S., Basham, P.W. (Eds.), Earthquakes at North Atlantic passive margins: Neotectonics and post-glacial rebound, NATO ASI Series. Springer, pp. 114–174.

Muller, K., Winsemann, J., Tanner, D.C., Lege, T., Spies, T., Brandes, C., 2021. Glacially Induced Faults in Germany. In: Steffen, H., Olesen, O., Sutinen, R. (Eds.), Glacially-Triggered Faulting. Cambridge University Press.

Nocquet, J.-M., 2012. Present-day kinematics of the Mediterranean: A comprehensive overview of GPS results. Tectonophysics 579, 220–242. https://doi.org/10.1016/j.tecto.2012.03.037.

Nocquet, J.-M., Calais, E., Parsons, B., 2005. Geodetic constraints on glacial isostatic adjustment in Europe. Geophys. Res. Lett. 32.

Ojala, A.E.K., Mattila, J., Hämäläinen, J., Sutinen, R., 2019. Lake sediment evidence of paleoseismicity: Timing and spatial occurrence of late- and postglacial earthquakes in Finland. Tectonophysics 771. https://doi.org/10.1016/j.tecto.2019.228227.

Peltier, W.R., 2004. Global Holocene Isostasy and the Surface of the Ice-Age Earth: The ICE-5G (VM2) Model and GRACE. Ann. Rev. Earth Planet. Sci. 32, 111–149.

Peltier, W.R., Drummond, R., 2008. Rheological stratification of the lithosphere: A direct inference based upon the geodynamically observed pattern of the glacial isostatic adjustment of the North American continent. Geophys. Res. Lett. 35 https://doi.org/10.1029/2008GL034586.

Pina-Valdés, J., Socquet, A., Beauval, C., Doin, M.-P., D'Agostino, N., Shen, Z.-K., 2022. 3D GNSS Velocity Field Sheds Light on the Deformation Mechanisms in Europe: Effects of the Vertical Crustal Motion on the Distribution of Seismicity. J. Geophys. Res. 127 https://doi.org/10.1029/2021JB023451.

Rollins, C., Avouac, J.-P., 2019. A Geodesy- and Seismicity-Based Local Earthquake Likelihood Model for Central Los Angeles. Geophys. Res. Lett. 46, 3153–3162. https://doi.org/10.1029/2018GL080868.

Rollins, C., Freymuller, J.T., Sauber, J.M., 2020. Stress Promotion of the 1958 M<sub>w</sub> ~ Fairweather Fault Earthquake ad Others in Southeast Alaska by Glacial isostatic Adjustment and Inter-earthquake Stress Transfer. J. Geophys. Res. 126 https://doi.org/10.1029/2020JB020411.

Ryan, W.B.F., Major, C.O., Lericolais, G., Goldstein, S.L., 2003. Catastrophic flooding of the Black Sea. Ann. Rev. Earth Planet. Sci. 31, 525–554. https://doi.org/10.1146/annurev.earth.31.100901.141249.

Sauber, J., Ruppert, N.A., 2008. Rapid Ice Mass Loss: Does It Have an Influence on Earthquake Occurrence in Southern Alaska? In: Freymuller, J.T., Haeussler, P.J., Wesson, R.L., Ekström, G. (Eds.), Active Tectonics and Seismic Potential of Alaska, Geophysical Monograph, vol. 179. American Geophysical Union. https://doi.org/10.1029/179GM21.

Sauber, J., Rollins, C., Freymuller, J., Ruppert, N., 2021. Glacially induced faulting in Alaska. In: Steffen, H., Olesen, O., Sutinen, R. (Eds.), Glacially-Triggered Faulting. Cambridge University Press.

Sauber, J.M., Molnia, B.F., 2004. Glacier ice mass fluctuations and fault instability in tectonically active southern alaska. Global Planet. Change 42, 279–293. https://doi.org/10.1016/j.gloplacha.2003.11.012.

Sella, G.F., Stein, S., Dixon, T.H., Craymer, M., James, T.S., Mazzotti, S., Dokka, R.K., Jan 2007. Observation of glacial isostatic adjustment in “stable” North America with GPS. Geophys. Res. Lett. 34 (2), L02306.

Steffen, H., Steffen, H., Tarasov, L., 2019. Modelling of glacially-induced stress changes in Latvia, Lithuania and the Kaliningrad District of Russia. Baltica 32, 78–90. https://doi.org/10.5200/baltica.2019.1.7.

Steffen, H., Olesen, O., Sutinen, R., 2021. Glacially triggering fault: a historical overview and recent developments. In: Steffen, H., Olesen, O., Sutinen, R. (Eds.), Glacially-Triggered Faulting. Cambridge University Press.

Stein, S., Liu, M., 2009. Long aftershock sequences within continents and implications for earthquake hazard assessment. Nature 462, 87–89. https://doi.org/10.1038/nature08502.
Stein, Seth, Geller, Robert J., Liu, Mian, 2012. Why earthquake hazard maps often fail and what to do about it. Tectonophysics 562–563, 1–25.

Sternai, P., Sue, C., Husson, L., Serpelloni, E., Becker, T.W., Willett, S.D., Faccenna, C., Di Giulio, A., Spada, G., Jolivet, L., Valla, P., Petit, C., Nocquet, J.-M., Walpersdorf, A., Castelltort, S., 2019. Present-day uplift of the European Alps: Evaluating mechanisms and models of their relative contributions. Earth Sci. Rev. 190, 589–604. https://doi.org/10.1016/j.earscirev.2019.01.005.

Townend, J., Zoback, M.D., 2000. How faulting keeps the crust strong. Geology 28 (5), 399–402.

Tuttle, M.P., Schwieg III, E.S., Campbell, J., Thomas, P.N., Sims, J.D., Lafferty III, R.H., 2005. Evidence for New Madrid Earthquakes in A.D. 300 and 2350 B.C. Seismol. Res. Lett. 76, 489–502.

van Balen, R.T., Bakker, M.A.J., Kasse, C., Wallinger, J., Woolderink, H.A.G., 2019. A Late Glacial surface rupturing earthquake at the Peel Boundary fault zone, Roer Valley Rift System the Netherlands. Quatern. Sci. Rev. 218, 254–266. https://doi.org/10.1016/s00531-020-01874-0.

Vanderberghe, D., Vanneste, K., Verbeeck, K., Paulissen, E., Buylaert, J.-P., De Corte, F., Van den Haute, P., 2009. Late Weichselian and Holocene earthquake events along the Geleen fault in NE Belgium: OSL age constraints. Quatern. Int. 199, 56–74. https://doi.org/10.1016/j.quaint.2007.11.017.

Vanneste, K., Verbeeck, K., Camlebeeck, T., Paulissen, E., Meghraoui, M., Renardy, F., Jongmans, D., Frechen, M., 2001. Surface-rupturing history of the Bree fault scarp, Roer Valley graben: Evidence for six events since the late Pleistocene. J. Seismolog. 5, 329–359.

Vanneste, K., Camelbeeck, T., Verbeeck, K., 2013. A Model of Composite Seismic Sources for the Lower Rhine Graben, Northwest Europe. Bull. Seismol. Soc. Am. 103, 984–1007. https://doi.org/10.1785/0120120037.

Štěpančíková, P., Hók, J., Nývlt, D., Dohnal, J., Sýkorová, I., Stemberk, J., 2012. Active tectonics research using trenching technique on the south-eastern section of the Sudetic Marginal Fault (NE Bohemian Massif, central Europe). Tectonophysics 485, 269–282. https://doi.org/10.1016/j.tecto.2010.01.004.

Štěpančíková, P., Fischer, T., Stemberk, J., Nováková, L., Hartvich, F., Figueiredo, P.M., 2019. Active tectonics on the Cheb Basin: youngest documented Holocene surface faulting in Central Europe. Geomorphology 327, 472–488. https://doi.org/10.1016/j geomorph.2018.11.007.

Štěpančíková, P., Rockwell, T.K., Stemberk, J., Rhodes, E.J., Hartvich, F., Luttrell, K., Myers, M., Táborík, P., Rood, D.H., Wechsler, N., Nývlt, D., Ortuno, M., Hók, J., 2022. Acceleration of Late Pleistocene activity of a Central European fault driven by ice loading. Earth Planet. Sci. Lett. 591 https://doi.org/10.1016/j. epsl.2022.117596.

Whitehouse, P., Gomez, N., King, M.A., Wiens, D.A., 2019. Solid Earth change and the evolution of the Antarctica Ice Sheet. Nat. Commun. 10 https://doi.org/10.1038/s41467-018-08068-y.

Wu, P., Johnston, P., 2000. Can deglaciation trigger earthquake in N. America? Geophys. Res. Lett. 27, 1323–1326.

Wu, P., Johnston, P., Lambeck, K., 1999. Post-glacial rebound and fault instability in Fennoscandia. Geophys. J. Int. 139, 657–670. https://doi.org/10.1046/j.1365-246x.1999.00963.x.

Zhao, S., Lambeck, K., Lidberg, M., 2012. Lithosphere thickness and mantle viscosity inverted from GPS-derived deformation rates in Fennoscandia. Geophys. J. Int. 190 https://doi.org/10.1111/j.1365-246X.2012.05454.x.

Zoback, M.D., Healy, J.H., 1992. In situ stress measurements to 3.5 km depth in the Cajon Pass Scientific Research Borehole: Implications for the mechanics of crustal faulting. J. Geophys. Res. 97, 5039–5057.