Numerical models including one or more faults in a rheologically stratified lithosphere show that climate-induced variations in ice and water volumes on Earth’s surface considerably affect the slip evolution of both thrust and normal faults. In general, the slip rate and hence the seismicity of a fault decreases during loading and increases during unloading. Here, we present several case studies to show that a postglacial slip rate increase occurred on faults worldwide in regions where ice caps and lakes decayed at the end of the last glaciation. Of note is that the postglacial amplification of seismicity was not restricted to the areas beneath the large Laurentide and Fennoscandian ice sheets but also occurred in regions affected by smaller ice caps or lakes, e.g. the Basin-and-Range Province. Our results do not only have important consequences for the interpretation of palaeoseismological records from faults in these regions but also for the evaluation of the future seismicity in regions currently affected by deglaciation like Greenland and Antarctica: shrinkage of the modern ice sheets owing to global warming may ultimately lead to an increase in earthquake frequency in these regions.
During the past two decades, much effort has been put into deciphering the interaction between climate changes and the deformation of the Earth’s lithosphere using numerical modelling. The starting point of most models was that topography created by tectonic processes will ultimately be the target of erosion and sedimentation processes, which will lead to the re-distribution of mass at the Earth’s surface. This mass re-distribution may, in turn, influence the localization and style of crustal deformation, as demonstrated by numerical modelling on different scales (e.g. Koons 1989; Beaumont et al. 1992; Avouac & Burov 1996; Batt & Braun 1999). Commonly, however, it has proven difficult to test the model predictions by field data, because cause and potential effect may be both temporally and spatially separated—in particular, on the scale of an orogen—and hence not easily identified.
On timescales of 104–105 years, a clear link between the slip history of active faults and climate-driven mass fluctuations on Earth’s surface was recently established by comparing predictions from numerical models with geological and palaeoseismological data (Hetzel & Hampel 2005; Hampel & Hetzel 2006; Hampel et al. 2007, 2009; Turpeinen et al. 2008). Such mass fluctuations include variations in the volumes of ice caps, glaciers or lakes and may considerably affect the slip behaviour of faults. The numerical models, which explicitly include one or several fault planes embedded in a rheologically layered lithosphere, showed that loading and flexure of the lithosphere by ice or water decreases the slip rates of nearby faults whereas unloading and rebound of the lithosphere accelerates the slip accumulation.
After a summary of the general setup and results of these previous models, this paper presents an overview of prominent examples from Scandinavia and the Basin-and-Range Province, for which numerical models showed that the postglacial increase in seismicity documented by field data can be explained by rapid changes in the volume of nearby ice caps or lakes. Finally, we describe regions where similar processes could have played a role and discuss the implications for regions that currently experience ice loss or lake regression.
2. General model setup and results
The principle setup of both the two- and three-dimensional numerical models includes a lithosphere that is divided into an elastic upper crust, a viscoelastic lower crust and a viscoelastic lithospheric mantle (figure 1). In the upper crust, one or more faults may be embedded, which may have variable dip and—in three-dimensional models—also variable strike. Fault slip is governed by a Mohr–Coulomb criterion. To simulate a contractional or extensional tectonic setting and to initiate slip on the fault, the model is shortened or extended by applying a velocity boundary condition (typically a few millimetres per year) at the model sides. If the models are applied to specific case studies, the shortening/extension rate is adjusted to velocities derived from space-geodetic measurements. The load, i.e. the ice or water volume, is represented by a pressure, which is applied to the top of the model and may vary both in space and time. Detailed information on the setup and the rheological parameters can be found in Hampel & Hetzel (2006) and Turpeinen et al. (2008).
The models are calculated using the commercial finite-element software ABAQUS. Each experiment starts with the establishment of isostatic equilibrium and an initial phase, during which the model is shortened or extended for usually several hundred thousand years. During this phase, the fault reaches a steady-state slip rate (figure 2a). Afterwards the cycle of loading and subsequent unloading begins, with the shortening/extension rate unchanged. In most experiments of our sensitivity analyses (Hampel & Hetzel 2006; Turpeinen et al. 2008; Karow 2009), the fault responds to loading and flexure of the lithosphere by a decreasing slip rate, whereas subsequent unloading and rebound accelerates fault slip (figure 2a). Parameters that have a major influence on the slip rate variations include the characteristics of the load (volume and spatial distribution), the viscosity of the lithospheric layers and the shortening/extension rate, which is applied as a boundary condition. In particular, the slip rate variations are considerably amplified with increasing magnitude of the load. Parameters of minor influence include the thickness of the lithosphere and the viscosity of the asthenosphere.
For both thrust and normal faults, the slip rate variations are caused by temporal changes in the differential stress, i.e. the difference between the maximum principal stress σ1 and the minimum principal stress σ3. For a fault located beneath the load, the differential stress (σ1−σ3) decreases during loading and flexure of the lithosphere and increases during unloading and rebound (figure 2b). As a consequence, the fault slip rate first decreases and subsequently increases. Remarkably, both thrust and normal faults show a similar slip pattern during the loading–unloading cycle even though the orientation of the principal stresses is different for contractional and extensional tectonic regimes (i.e. σ1 is horizontal in contractional settings, whereas σ1 is vertical in extensional regimes). The underlying reason is that—for both settings—the horizontal stress changes, which are caused by flexure and rebound of the lithosphere, exceed the vertical stress changes, which are caused by addition and removal of the load. Hence the differential stress shows a similar temporal evolution for both thrust and normal faults.
3. Case studies
One region where a close spatio-temporal correlation between melting of a large ice sheet and the occurrence of palaeo-earthquakes was already recognized in the 1970s is the Lapland Fault Province in northern Scandinavia (figure 3a; Lundquist & Lagerbäck 1976; Mörner 1978; Lagerbäck 1979). In this region—commonly considered as a tectonically quiet continental shield—up to 15 m high scarps along faults up to 150 km long document the occurrence of palaeo-earthquakes with magnitudes as large as Mw≈8 (e.g. Arvidsson 1996; Dehls et al. 2000; Lagerbäck & Sundh 2008). Most of the faults have a reverse sense of slip and strike NNE–SSW. Prominent examples include the Pärvie fault in northern Sweden (figure 3b,c) and the Stuoragurra fault in Norway. As the fault scarps offset deposits from the last glacial period, they must be of postglacial age. Trenching and dating of earthquake-induced soil liquefaction features and landslides revealed that the palaeo-earthquakes in the Lapland Fault Province clustered at approximately 9000 yr BP (Lagerbäck 1992; Mörner 2005), which coincides in space and time with the deglaciation of the region (e.g. Lundquist 1986).
To investigate whether the palaeo-earthquakes in Scandinavia were triggered by the melting of the more than 2 km thick Fennoscandian ice sheet, we used a two-dimensional model, in which the width of the ice sheet, its temporal evolution as well as the tectonic shortening rate were adjusted to northern Scandinavia (for details see Turpeinen et al. 2008). The model results showed that faulting is suppressed during the entire lifetime of the Fennoscandian ice sheet, regardless of temporal variations in its thickness between 700 and 3400 m (figure 4). Melting of the ice sheet triggers a strong slip rate increase: between approximately 10.5 and 8 kyr, the fault accumulates approximately 15 m of slip. This is in excellent agreement with the timing and the displacement of the palaeo-earthquakes in the Lapland Fault Province. The model results further reveal that most of the slip is accumulated on the fault before deglaciation is completed, which is in accordance with the palaeoseismological record. Note that the model does not predict the number of earthquakes by which the total displacement is achieved. In the Lapland Fault Province, some large faults like the Pärvie and Lansjärv faults probably ruptured in a single event (Arvidsson 1996; Dehls et al. 2000). In summary, our models support the view that deglaciation caused the increased seismicity in northern Scandinavia. The same may be true in central and southern Scandinavia, where palaeo-earthquakes also coincided in space and time with the deglaciation. In central Sweden an increase in seismicity is documented at approximately 10 500 varve yr BP (Mörner 2005), whereas palaeoseismological investigations in southern Sweden revealed a cluster of earthquakes between approximately 12 000 and 12 400 varve yr BP. Hence, it appears that deglaciation-induced seismicity migrated northward together with the margin of the melting Fennoscandian ice sheet.
(b) Basin-and-Range Province
In contrast to Scandinavia, the Basin-and-Range Province in the western United States is characterized by active east–west extension at a rate of several millimetres per year (e.g. Wernicke et al. 2000), which led to the formation of alternating normal-fault bounded ranges and sedimentary basins. For several of the many normal faults in the Basin-and-Range Province, detailed geological and palaeoseismological records constrain the slip histories on time scales of up to 105 years. Interestingly, some of these records revealed considerable variations in the fault slip rates on a 104 years time scale (e.g. Wallace 1987; Byrd et al. 1994; Friedrich et al. 2003; Wesnousky & Willoughby 2003; Wesnousky et al. 2005). Such highly variable faulting rates are difficult to reconcile with quasi-periodic earthquake recurrence models and lead to a new recurrence model that combines earthquake clustering on short time scales with uniform long-term strain accumulation (Wallace 1987). It is intriguing, however, that several faults including the prominent Wasatch and Teton faults show an increase in their slip rate at the end of the last glacial period (Byrd et al. 1994; McCalpin 2002; Ruleman & Lageson 2002). Although the Basin-and-Range Province was not covered by a large ice sheet like Scandinavia, it experienced considerable loading by two large lakes—Lake Bonneville and Lake Lahontan—and an ice cap centred on the Yellowstone Plateau during the last glacial period. However, the potential relationship between the glacial loads and accelerated fault slip is not as obvious as in Scandinavia owing to the active crustal extension. In the following, we present two case studies, for which well-documented fault slip variations could be linked to changing ice and water volumes (Hetzel & Hampel 2005; Hampel et al. 2007; Karow & Hampel in press).
The first case study is dedicated to the approximately 350 km long Wasatch fault in the eastern Basin-and-Range Province and its response to the regression of Lake Bonneville (figure 5; Hetzel & Hampel 2005; Karow & Hampel in press). This predecessor of the Great Salt Lake occupied an area of approximately 52 000 km2 and reached a maximum depth of 350 m. Well-preserved shorelines on islands and the flanks of mountain ranges (figure 5b,c) reveal three highstands (Gilbert 1890) at 17 500 calibrated 14C yr BP (17.5 kyr; Bonneville highstand), 16.7 kyr (Provo highstand) and 11.4 kyr (Gilbert highstand), respectively (Oviatt et al. 1992). As first noted by Gilbert (1890), the shorelines in the former centre of Lake Bonneville occur at an elevation 60–70 m higher than the shorelines along the lake periphery, which documents that the lake regression caused pronounced isostatic uplift (Bills et al. 1994). In addition to Lake Bonneville, the area east of the Wasatch fault was loaded by glaciers in the Uinta and Wasatch Mountains (figure 5d,e), which locally reached a thickness of several hundred metres (Atwood 1909; Madsen & Currey 1979).
A wealth of palaeoseismological data exists for the different segments of the Wasatch fault and the data sets from the Salt Lake City and Brigham City segments even extend beyond the Holocene (e.g. Personius 1991; McCalpin et al. 1994; McCalpin & Nishenko 1996; McCalpin 2002; McCalpin & Forman 2002). The Salt Lake City segment experienced a phase of seismic quiescence lasting from at least 26 kyr to the end of the Bonneville highstand, which was followed by seven earthquakes. The first of these earthquakes occurred at approximately 17 kyr, i.e. shortly after the rapid regression of the lake; the six other seismic events followed between approximately 9 kyr and present. Similar to the Salt Lake City segment, the Brigham City segment was ruptured by an earthquake that occurred after the Bonneville flood but before the end of the Provo highstand. Another palaeo-earthquake probably occurred between 9 and 13 kyr BP, followed by a sequence of five earthquakes during the Holocene. From the cumulative post-Bonneville displacements on the Salt Lake City and Brigham City segments an average slip rate of approximately 1 mm yr−1 since regression of Lake Bonneville was inferred (e.g. Friedrich et al. 2003). In contrast, the average slip rate of the central Wasatch fault during the last 130 kyr was only 0.1–0.6 mm yr−1, as derived from tectonically offset Bull Lake glacial deposits and pre-Bonneville lacustrine deposits. This implies that the post-Bonneville slip rate is about twice as high as the average slip rate on a 105 years time scale. Slip accumulation on the Wasatch fault varied also along-strike because the cumulative displacement since the regression of Lake Bonneville is 18–20 m on the Brigham City segment but only 5–10 m on the southern Nephi segment.
Using two-dimensional models including a normal fault, Hetzel & Hampel (2005) showed that the rebound induced by regression of the lake and melting of glaciers increases the slip rate of the model fault by temporally altering the stress field of the crust. To investigate potential along-strike slip variations that may result from the spatial distribution of Lake Bonneville and the valley glaciers in the Uinta and Wasatch mountains, Karow & Hampel (in press) recently created a three-dimensional numerical model, which takes the spatial distribution of the lake and glaciers into account. This model shows a maximum rebound-induced uplift of 69 m, which agrees well with the uplift derived from deformed lake shorelines, and reveals substantial along-strike variations in fault slip (figure 6). During rise of the lake level, from 35 to 17.5 kyr, most slip is accumulated in the southern part of the model fault. From 17.5 kyr to present, the northern and central segments experience the most slip (figure 6). In particular, the model predicts a pronounced slip rate increase on the northern and central segments shortly after the Bonneville flood (figure 6b), which agrees well with the occurrence of the first post-Bonneville earthquakes on the Brigham City and Salt Lake City segments approximately 17 kyr ago (McCalpin 2002; McCalpin & Forman 2002). In accordance with the geological and palaeoseismological data (McCalpin & Nishenko 1996; McCalpin 2002; Friedrich et al. 2003), the slip rate of the model fault since 17.5 kyr is higher than its long-term slip rate. Finally, the model shows that the central Brigham City and Salt Lake City segments accumulate about three times more displacement than the southern Nephi segment, which agrees with the observed along-strike variations in fault slip.
The second case study from the Basin-and-Range Province sheds light on the influence of the Yellowstone ice cap on the Teton normal fault (figure 7). The Teton Range in the fault footwall is famous for its spectacular glacial geomorphology and high relief. Its highest peak, the Grand Teton, stands approximately 2400 m above Jackson Hole, the hanging wall of the fault (figure 7b,c). Both elevation and relief along the range decrease along-strike towards the north and south (figure 7b). Maximum estimates for the slip rate of the Teton fault on a 106 years time scale range between 0.5 and 1.2 mm yr−1 (Love 1977; Smith et al. 1993; Byrd et al. 1994). Fault scarps, which offset deposits from the last glacial period (Pinedale), provide evidence for its seismic activity during the Late Pleistocene and Holocene (figure 7d). Vertical offsets derived from these scarps reach approximately 30 m in the central part of the fault and diminish to the north and south, thus mimicking the along-strike variations in elevation and relief of the range. Assuming that the scarps formed since approximately 15 kyr, a vertical slip rate of approximately 2 mm yr−1 was derived for the central Teton fault (Machette et al. 2001). Trenching across the southern Teton fault revealed two Holocene earthquakes at approximately 8 and 5 kyr with a total vertical offset of approximately 4 m (Smith et al. 1993; Byrd et al. 1994). As the surface offset near the trench site is approximately 15 m, approximately 11 m of vertical displacement must have been produced between 8 kyr and the onset of deglaciation (Lageson et al. 1999). Hence approximately 70 per cent of the postglacial slip on the Teton fault occurred before 8 kyr but after the end of the Pinedale glaciation (16–14 kyr; Licciardi et al. 2001; Licciardi & Pierce 2008). In the valleys adjacent to the Grand Teton and Mt Moran (figure 7b,c), the Pinedale glaciers attained a maximum thickness of approximately 1000 m (Love et al. 2003). Additional loading of the Teton region was caused by the Yellowstone ice cap, which covered an area of approximately 16 500 km2 and reached a surface altitude of approximately 3400 m (Love et al. 2003). Terminal moraines indicate that their southernmost lobes reached into northern Jackson Hole, where the melting of the ice formed Jackson Lake (figure 7c).
To evaluate whether the melting of the Yellowstone ice cap and the Teton valley glaciers was able to trigger the postglacial slip rate increase on the Teton fault, a three-dimensional numerical model is required because the ice cap is mostly located north of the Teton fault (figure 7a). The results of our three-dimensional model (Hampel et al. 2007) show that the ice load caused more than 80 m of subsidence beneath the ice cap centre 22 kyr ago (figure 8a). Growth and subsequent melting of the ice causes considerable slip rate variations on the model fault. The slip rate at the fault centre on the model surface decreases markedly during growth of the ice cap. During deglaciation, the slip rate increases by a factor of approximately 5 with respect to the steady-state rate and reaches a value of 2.33 mm yr−1 (figure 8b). The model further reveals pronounced variations in the postglacial slip along-strike of the fault, because the influence of the Yellowstone ice cap on fault slip increases from south to north (cf. Hampel et al. 2007). To compare our results with the palaeoseismological data, we calculated the ratio between the slip accumulated from 16 to 8 kyr and the slip accumulated from 8 kyr to present. For the model without loading, the slip accumulated in both time intervals is equal. In contrast, loading by the Teton valley glaciers and the Yellowstone ice cap leads to a slip distribution with two-thirds of the slip occurring between 16 and 8 kyr (figure 8b), which agrees with the palaeoseismological data. Our model further predicts that a considerable amount of slip occurred between 16 and 14 kyr. Unfortunately, the existing palaeoseismological data do not allow a detailed comparison with this model prediction because the timing of palaeo-earthquakes between 16 and 8 kyr is still unknown.
4. Implications for other formerly glaciated mountain ranges and for regions currently experiencing ice mass loss
Spatial and temporal slip rate variations are often taken as indicative of processes during fault growth, interaction among faults and fault linkage (Anders & Schlische 1994; Cowie 1998; Roberts & Michetti 2004). In particular, postseismic relaxation and reloading may affect the future seismic activity of fault networks (Kenner & Simons 2005; Meade & Hager 2004). Temporal slip rate variations may further arise from changes in the friction coefficient of the fault, which may in turn be caused by climate-controlled changes in pore pressure (Chéry & Vernant 2006). In some cases, an apparent increase in seismicity during the Holocene may simply be the result of a sampling bias towards faults with higher Holocene slip rates because these faults have a clearer morphological expression than faults with low slip rates (Nicol et al. 2009). However, if several faults show a synchronous increase in their slip rate with a timing correlated with major climate transitions, an external driving force from changing loads on Earth’s surface should be considered. During the Last Glacial Maximum, numerous mountain ranges worldwide were covered by ice caps or valley glaciers, for example, the New Zealand Alps, the North American Cordillera, the southern Andes and the European Alps (e.g. Ehlers & Gibbard 2004). It can thus be expected that more faults than the examples presented above were affected by postglacial unloading and rebound. For the European Alps, which were covered by up to 2 km of ice during the Last Glacial Maximum, an integrated field and modelling study shed light on the formation of unusual uphill-facing fault scarps occurring on hillslopes in the central Swiss Alps (Ustaszewski et al. 2008). These up to 10 m high scarps run parallel to the axes of formerly glaciated valleys and locally offset deposits from the last glaciation. Finite-element experiments showed that their formation is linked to the deglaciation of the valleys, which caused uplift of the valley floors with respect to the ridge crests. A recent three-dimensional numerical models suggests that deglaciation of the European Alps after the Last Glacial Maximum caused more than 130 m of rebound, which still continues today at a rate of approximately 0.3 mm yr−1 (Norton & Hampel submitted). Based on these model results, we speculate that the increase in seismicity after the Last Glacial Maximum inferred from palaeoseismological records in the French and Swiss Alps (Beck et al. 1996; Becker et al. 2005 and references therein) may be triggered by postglacial unloading and rebound. Another example from this region is a normal fault in the southern Upper Rhine Graben, which was related to the Basel earthquake in 1356. Palaeoseismological investigations suggest that five earthquakes with a cumulative slip of approximately 3 m occurred on this fault during the last 13 kyr, which is equivalent to a vertical displacement rate of 0.27 mm yr−1 (Ferry et al. 2005). This rate is by a factor of three higher than the average faulting rate of 0.11 mm yr−1 over the last 0.7 Myr. The increased postglacial slip accumulation of the fault may be related to both the deglaciation of the Alps and of glaciers in the southern Black Forest, Germany, and the southern Vosges, France (cf. Ehlers & Gibbard 2004 and references therein).
With respect to currently glaciated regions, our models indicate that the currently low level of seismicity in Greenland and Antarctica is caused by the presence of the large ice sheets. They further imply that the frequency of earthquakes should increase in the future if the ice continues to melt. This effect may be important even on time scales of 10–100 years, as shown by studies in south-central Alaska where the level of seismicity during the last decades was modulated by glacial mass fluctuations (e.g. Sauber & Ruppert 2008). In particular, the number and size of earthquakes increased during a large-scale redistribution of ice, which occurred between 1993 and 1995 and unloaded an area near the Bering glacier in the eastern Chugach Mountains.
Numerical models show that mass fluctuations on Earth’s surface affect the slip evolution of both normal and reverse faults and provide quantitative predictions on this type of interaction between tectonics and climate. In general, a low level of seismicity can be expected in the presence of a load whereas slip is accelerated during unloading and rebound. Hence glacial loading and subsequent postglacial unloading provide a feasible mechanism to explain the synchronous increase in seismicity reported from formerly glaciated regions worldwide. Parameters exerting a major control on fault behaviour include the characteristics of the load (volume, spatial distribution), the viscosity of the lower crust and the lithospheric mantle and the rate of shortening or extension. Based on the model results, it is expected that the decay of the modern ice sheets in Greenland and Antarctica will ultimately increase the level of seismicity in these regions.
Over the last years, a number of scientific partners have contributed to the topic presented in this overview article, either by contributions to previous publications or inspiring discussions. In particular, we thank A. Densmore, M. Ustaszewski, A. Pfiffner, G. Roberts, J. Sauber and I. Stewart as well as our former MSc student H. Turpeinen and doctoral student T. Karow. R. Lagerbäck kindly provided the two pictures from the Pärvie fault. We thank an anonymous reviewer and the editor B. McGuire for their positive reviews. Funding by the German Research Foundation (DFG) within the framework of an Emmy-Noether fellowship to A.H. (grant HA 3473/2-1) is gratefully acknowledged.
One contribution of 15 to a Theme Issue ‘Climate forcing of geological and geomorphological hazards’.
- © 2010 The Royal Society