Earthquake-induced transformation of the lower crust.

The structural and metamorphic evolution of the lower crust has direct effects on the lithospheric response to plate tectonic processes involved in orogeny, including subsidence of sedimentary basins, stability of deep mountain roots and extension of high-topography regions. Recent research shows that before orogeny most of the lower crust is dry, impermeable and mechanically strong 1 . During an orogenic event, the evolution of the lower crust is controlled by infiltration of fluids along localized shear or fracture zones. In the Bergen Arcs of Western Norway, shear zones initiate as faults generated by lower-crustal earthquakes. Seismic slip in the dry lower crust requires stresses at a level that can only be sustained over short timescales or local weakening mechanisms. However, normal earthquake activity in the seismogenic zone produces stress pulses that drive aftershocks in the lower crust 2 . Here we show that the volume of lower crust affected by such aftershocks is substantial and that fluid-driven associated metamorphic and structural transformations of the lower crust follow these earthquakes. This provides a 'top-down' effect on crustal geodynamics and connects processes operating at very different timescales.

The structural and metamorphic evolutions of the lower crust are key elements in the dynamics of the lithosphere. Frequent observations of fluid-induced metamorphism associated with ductile deformation along shear zones on scales ranging from millimeters to kilometers inspired early models of the lithosphere such as the ‘jelly-sandwich’ model3. In this model, the lower crust is assumed to be wet and mechanically weak, and plate tectonic stress is transmitted through the brittle upper crust and a strong upper mantle. This model was challenged4 with the argument that a strong lower crust is essential for the survival of thick mountain roots and high mountains. The lower crust is dominated by granulite facies rocks of mafic to intermediate composition5 and such rocks will be nominally dry at normal steady state geothermal gradients for a wide range of crustal heat flow and heat production conditions1. Hence, the rheology of the lower crust prior to an orogeny will in most cases be controlled by the properties of dry mineral assemblages dominated by plagioclase, pyroxene, garnet and olivine, with plagioclase being the most abundant phase. This is consistent with the estimated viscosity of the lower crust (≥ 1024 Pa·s) required to generate the crustal support needed for intraplate seismicity such as the 2001 MW 7.6 Bhuj earthquake in western India6. Recent modeling7 furthermore suggests that the Indian lower crust remains strong beneath the entire southern half of the Tibetan plateau.

Observations on structural and metamorphic transformation of initially dry lower crust during orogenic events indicate an early stage involving seismic failures8–10. Metamorphism and shear zone development then follow in the wake of lower crustal earthquakes. These observations raise an enigma that has so far been unresolved, because frictional failure of dry rocks at the confining pressures of the lower crust (≥ 1 GPa) requires differential stress levels exceeding 2 GPa11. Although dry plagioclase-dominated rocks deforming by dislocation creep can in theory develop extremely high differential stress at lower crustal temperatures and high strain rates, the stress level that can be sustained over orogenic time scales for reasonable strain rates in coherent crustal volumes will be far below what is required for brittle faulting (< 1 GPa, see Methods). Deep crustal earthquakes occurring under constant loading in intact rocks thus seem to require a local weakening mechanism.

During subduction of the Indian plate under south Tibet, earthquakes occur at 60-100 km depth (Fig.1), but are confined to regions very close to the Mohorovičić discontinuity (Moho) at temperatures below 600°C7. Interestingly, the pressure and temperature conditions in the region where these earthquakes nucleate overlap the conditions at which serpentine breaks down to produce hydrous fluids in mantle rocks. Fluid production near the Moho can both reduce effective pressures and weaken the crust and mantle by mineral transformation processes and thus be a plausible explanation for the observed seismic activity.

Figure 1

Earthquakes and aftershocks in the lower crust.

a, Schematic representation of earthquakes and aftershocks for the India-Tibet continent-continent collision. b, A generic subduction plate boundary geometry. Each major earthquake generates a cloud of aftershocks, some of which are in the lower crust (many aftershocks in the upper crust are not shown in these simple diagrams). These aftershocks create pathways for fluids (blue arrows in insets), allowing partial hydration and metamorphism of the strong and dry granulites into wet and weaker eclogites and amphibolites. This process also facilitates the development of shear zones in the continental lower crust. For the subduction geometry, fluids could originate from the slab below or from the upper crust above. For the continental collision below the Himalayas, fluids introduced to the subducted Indian plate could originate from the dehydration of serpentine rocks below. The 2001 Bhuj6 and 2015 Gorkha12 earthquakes both have aftershock ‘clouds’ propagating down to the lower crust, potentially allowing downward migration of fluids from the upper crust.

In the absence of such local weakening mechanisms, seismic deformation in the lower crust may be driven by transient ‘stress pulses’ 6,10,13. Here we propose that regular earthquakes in the brittle upper crust provide a natural mechanism for sustained generation of stress pulses and associated seismic failures in the lower crust. During the occurrence of large earthquakes the strain rates around and below the rupture area increase by many orders of magnitudes. A representative strain accumulation of 15 mm per year across a width comparable to a geodetic locking depth of 15 km corresponds to an interseismic strain rate of 3·10-14/s. In contrast, seismic slip velocities of 1-10 m/s across rupture localization width of 1-10 mm lead to seismic strain rates of 103-104/s. Such large co-seismic jumps can explain a transient increase in seismic rupture within the lower crust.

Observational evidence for very high, short-lived, stresses in the lower crust come from the occurrence of fossil earthquakes. A recent study of the Woodroffe Thrust located within the Musgrave Block in Central Australia14 documents the formation of large volumes of pseudotachylytes in completely dry lower crustal granulites. Stress levels exceeding 0.5 GPa have also been inferred from lower crustal earthquakes leading to pseudotachylyte formation in gabbros and ultramafic rocks in the Alpine subduction complex of Corsica, France15.

Simulated deformation on faults using various versions of rate- and state-dependent friction models show that large earthquake slip penetrates into the nominally stable deeper region16,17. Simulations of aftershocks in a viscoelastic damage model consisting of a brittle upper crust over a lower crust with power-law viscosity constrained by laboratory experiments indicate that the hypocenters of early aftershocks are significantly deeper than the regular seismogenic zone2. Depending on model parameters and thermal gradients, the maximum depth of the early aftershocks can approach twice that of the usual seismicity. Details of these results depend on the employed constitutive laws and parameters, but lower crustal aftershocks are generic outcomes of the high strain rates generated by large earthquakes at the bottom of the seismogenic zone.

Below we use basic scaling relations to demonstrate that observed worldwide earthquake activity in the regular seismogenic zone of subduction zones and seismically active continental regions is expected to produce considerable fracture area and rupture zone volume in the lower crust. This, in turn, generates transient pathways for fluids from the wet upper crust above, or the slab below, to the dry lower crust (Fig. 1). Fluids play a key role in the long-term evolution of the lower crust1. As an example, we describe earthquake-triggered eclogite-facies metamorphism and shear zone development of lower crustal granulites from the Bergen Arcs in Western Norway (Fig. 2A). The observations highlight the close association between earthquakes, fluid migration and transformation processes in the lower crust.

Figure 2

Fossil earthquakes in the Bergen Arcs.

a, Map of the Northwestern part of Holsnøy island in the Bergen Arcs, Western Norway (modified from Ref. 20) showing the location of pseudotachylytes (red stars) recording numerous fossil earthquakes near the transition between 940 Ma old dry lower crustal granulites and hydrated 430 Ma old eclogites. b, Offset of a pyroxene-rich seam by a single lower crustal earthquake. The slip surface is decorated with a melt layer, (pseudotachylyte) indicative of seismic slip. An offset of 1.7 m corresponds to a fossil earthquake with M ≥ 7 based on the scaling relations provided in Ref. 21. c, Higher resolution image of the centimeter-thick melt layer between the two red dashed lines. Red arrows indicate the sense of slip.

The Bergen Arcs represent a series of thrust sheets where granulite facies remnants of Proterozoic lower crust recrystallized to an anhydrous mineralogy at 930 Ma18. During the Caledonian continent-collision between Laurentia and Baltica between 420 and 440 Ma18, fluid-induced metamorphic transformations formed eclogites and amphibolites in shear zones, breccias and along fractures. The estimated eclogitization conditions are ca. 670-690°C and 2.1-2.2 GPa19. Pseudotachylytes, fine-grained or glassy fault rocks believed to reflect earthquake related frictional melting, are abundant on faults where granulites facies rocks experienced Caledonian retrograde metamorphism (Fig. 2A). Such faults show single rupture displacements reaching 1.7 m (Fig. 2B), corresponding to an earthquake exceeding magnitude 7. Single pseudotachylyte veins range in thickness from millimeters to a few centimeters (Fig. 2C), and also occur as a thin ‘matrix’ between rotated blocks of brecciated granulite that sometimes cover areas exceeding 100 m2. Microstructures developing in the fault wall rocks display intense fragmentation without significant shear strain22, followed by healing processes through grain growth and formation of eclogite facies minerals, including hydrous phases such as amphibole, mica and clinozoisite. Infiltration of hydrous fluids was thus directly associated with seismic slip.

A significant rheological weakening associated with formation of the fine-grained and hydrous eclogite often leads to development of ductile shear zones in areas initially deformed by brittle failure. Relict pseudotachylytes can occasionally be observed ‘floating’ in the shear zones, providing unambiguous evidence for ductile deformation predated by brittle failure of granulite facies rocks (Fig. 3). In the following, we explore the feasibility that lower crustal earthquakes, such as those described, are aftershocks triggered by stress pulses generated by mainshocks in the normal seismogenic regime of a plate boundary.

Figure 3

Transformation of the lower crust.

Offset of a pseudotachylyte by a shear zone where dry granulite rocks are transformed into wet eclogites. The earthquake occurred before the eclogitization and likely created the pathways for fluids that triggered rock transformation. Note the remains of the pseudotachylyte inside the eclogite shear zone.

Basic seismological scaling relations provide an order-of-magnitude estimate of the lower crustal rock volume affected by aftershocks. We demonstrate that this is significant, with conservative parameter values and ignoring probable contributions from penetration of large mainshocks into the lower crust as well as ductile/thermal instabilities16,17.

Lower crustal earthquakes are not expected to occur repeatedly in the same location because rock melting and subsequent solidification is a strengthening process23. This is consistent with observations of distributed “fields” of pseudotachylytes (Fig. 2A), each associated with a single event. The total volume of rock damage V produced by crustal earthquakes in the magnitude range M1 < M < M2 is then given by Here A(M) is the rupture area of an earthquake with magnitude M, t(M) is the damage zone thickness around the rupture area and n(M) is the event density given by the Gutenberg-Richter relationship where a and b are empirical constants. The scaling relations for A(M) and t(M) are found using basic theoretical relations from fracture mechanics24,25 and empirical relation between the magnitude and potency of earthquakes26,27 (see Methods). An explicit relationship for the total volume of damage produced by earthquakes in the considered magnitude range can be expressed as: where α and β are positive parameters that account for the combined scaling relations of A, t and n with M24–27. By using observationally-constrained parameters, the volume of rock damaged by crustal earthquakes is estimated to be V =1.2·10−5 km3 per year per km2 of the Earth’s surface in the seismically active region (see Methods). For the lower magnitude limit we use M1 = 0 with slip distance of the order of the grain size of granulites. For an upper limit relevant to the lower crust we take M2 = 8.3 since the largest subduction zones events can have M ≥ 9.5, and the largest aftershock magnitude is typically ~1.2 units below that of the mainshock28,29. The parameters of the Gutenberg-Richter relationship are taken from recent analysis of global earthquakes with depths less than 70 km30. Analysis of earthquake clusters indicates29,30 that ~75% of all events with M > 0 are aftershocks. We therefore use 75% of the observed intensity of events30 to estimate the average annual production of damaged rock volume by aftershocks with 0 < M < 8.3. Based on previously conducted simulations2 we assume conservatively that 1% of the aftershock population is in the lower crust. The estimated annual production of rupture zone volume in the lower crust is then 1.2∙10-7 km3/yr per km2 area at the Earth surface.

For a concrete example, Western Norway was an active subduction zone during the Caledonian orogeny for >2 ·106 yr. Based on the estimates above, the total seismically damaged volume in the lower crust of Western Norway is estimated to be 0.24 km3 per km2 area at the Earth surface. For a 20 km thick lower crust, this implies a rupture zone volume exceeding 1.2% of the total lower crustal volume. The Bergen Arcs example demonstrates that in the presence of fluids, lower crustal earthquakes initiate metamorphism of rock volumes typically 1-2 orders of magnitude larger than that of the rupture zone (i.e. a 0.1-1 m thick eclogite forming around 1 cm thick rupture zone). Hence the overall process can alter a large fraction of the lower crust. The sensitivity of the results to input parameters is discussed in the Methods section and shows that the 1.2% estimate of lower crust volume damaged by earthquakes is based on conservative values of input parameters and is likely to be higher.

Our results indicate that aftershocks triggered by major earthquakes in the regular seismogenic zone have the potential to initiate pervasive transformation of the lower crust on a timescale of 106 years. Direct recording of a transient deepening of early aftershocks requires a dense observational network around large mainshock ruptures. Although this situation is not often met, such lower crust aftershocks are sometimes observed. Recent examples of deep aftershocks include the 2001 M 7.6 Bhuj intraplate earthquake in India where aftershocks occurred down to Moho depths6, and the 2015 M 7.8 Gorkha earthquake in Nepal where the hypocenter occurred near the Main Himalayan Thrust at a depth of 10-15 km12 and aftershocks penetrated the Indian crust to a depth exceeding 30 km (Fig. 1A).

Earthquakes in dry lower crust driven by stress pulses generated in the seismogenic zone have a number of important consequences, many of which are illustrated by the Bergen Arcs example. The most important is arguably the associated increase in permeability which may connect the dry lower crust to an external fluid reservoir. In the Bergen Arcs, pseudotachylyte formation is always associated with influx of hydrous fluids21. Fluid driven metamorphic reactions are fast due to the metastable state of the granulite facies rocks1, leading to a profound reduction in rock strength and the development of shear zones and ductile deformation at lower stress levels31,32. The positive feedbacks between fluid introduction, weakening, and shear zone development eventually produce a complete transformation of large volumes of lower crust from a low density, dry and strong lithology, to a high density, wet and weak one. Thus even in situations where aftershocks directly affect only limited volumes of lower crust, they may start a series of fluid-induced transformation processes which can effect far bigger volumes.

The ‘top-down’ control on lower crustal evolution presented here challenges the traditional ‘bottom-up’ view where deep shear zones are assumed to control the spatial distribution of faults above the brittle-ductile transition. Generation of deep crustal shear zones as a response to weakening induced by pre-existing faults triggered by stress pulses generated by shallower earthquakes may also explain observed fluids with meteoric and other upper crustal signatures, such as the presence of hydrocarbons, in shear zones formed below the brittle-ductile transition33.

Methods

Controls on lower crustal stress levels

Flow law creep parameters for synthetic plagioclase aggregates for variable water contents were presented by Rybacki and Dresen34. Extended Data Fig.1 shows the relations between temperature, strain rate and differential stress for water poor plagioclase aggregates based on these experimental data. At 660-680°C, the estimated temperature during seismic faulting in the Bergen Arcs, the maximum differential stress developing for a strain rate of 10-14 s-1 would be ca. 0.3 GPa and clearly insufficient to cause brittle failure. Such strain rates would probably only apply within zones where strain was already localized. In a coherent crustal volume, the internal strain rate would be far less. Craig et al7., estimate strain rates at around 10-16 s-1 for the subducting Indian plate beneath eastern Tibet. Even if local strain rate increases should be able to push stress levels beyond what would be expected in coherent crustal slabs, or if the crustal temperature was lower than for the Bergen case so that the effective viscosity was higher, stresses much higher than 1 GPa would be unrealistic due to the onset of Peierls creep35, or even cataclastic flow36.

Extended Data Figure 1

Rheology of dry anorthite.

Shear stress versus temperature diagram contoured with respect to strain rate.

Volume of damage produced by earthquakes

The step-by-step derivation of Equation (3) and the calculation of the volume of rock damaged by earthquakes in the crust, V, are detailed here. The total volume of rock damage produced by crustal earthquakes in the magnitude range M1 < M < M2 is given by

Each term in the integral is described below.

The density of events n(M) in (S1) is provided by the Gutenberg-Richter relationship where a and b are empirical constants.

Assuming that each earthquake is approximately a circular rupture with radius r and surface area A = πr2 that sustains a uniform strain drop Δε, the seismic potency, P0, (moment/rigidity) is given by24,

The seismic potency and magnitude of earthquakes spanning a relatively small range (≤ 4) of magnitudes are related empirically by a relation of the type26,27, with potency units in km2·cm and constants d and e.

Combining Eqs. (S3) and (S4), the radius of an earthquake with magnitude M in unit of km is: where the 10−5 factor account for unit conversion of P0 from km2·cm to km3. Using (S5) for r and assuming again a circular crack such that A(M) = πr2, one obtains

The thickness of the damage zone, t, is expected from fracture mechanics to scale linearly with the rupture radius25. where the constant γ is proportional to the dynamic stress intensity factor and the ratio of stress drop over strength drop. For standard rupture velocity of 0.9 times the Rayleigh wave speed and a stress drop that is 10% of the strength drop, γ~10−2 − 10−3 (situations with relatively high initial stress leading to higher stress drop give higher γ values). Here we conservatively use γ = 0.002.

Combining (S2), (S6) and (S7) leads to

Integrating (S1) using (S8) gives an explicit relationship for the total volume of damage produced by earthquakes in the form where α = (d − b)ln(10) and are positive parameters and V is in unit of km3 per year per km2 of the Earth’s surface in a seismically active region.

Observed b-values are typically around 1 while a-values vary significantly with space and time. We focus on deformation in subduction zones and seismically active continental regions and use average a and b values based on analysis of global earthquakes with depth shallower than 70 km during 1975-2015 in the Northern California Earthquake Data Center (NCEDC) catalog30. From Fig. 1 of Ref. [30] and a b-value of 1, a representative intensity of events with M > 0 (10) in active subduction zones is 6/(yr·km2), corresponding to a = log10(6) = 0.78. This value is conservative since the NCEDC catalog does not include numerus small events buried in the coda of larger events and noise. The results of Ref. [30] indicate that about 50% of M > 4 events in subduction zones and continental transform regions are aftershocks. The fraction of events that are aftershocks increases as the event magnitude decreases29,30, so we assume that ~75% of all events with M > 0 are aftershocks. The event intensity 10 used below to estimate damaged rock volume by aftershocks is reduced accordingly by a factor of 0.75 from 6 to 4.5.

Analysis of earthquakes in southern California recorded by the regional network and borehole sensors indicates26 that d = 1 and e = −4.7 for M < 3.5 while d = 1.34 and e = −5.22 for M > 3.5. Similar constants characterize earthquakes in other locations27. Inserting into Eq. (S8) these constants d and e for 0 < M < 3.5 and 3.5 < M < 8.3 magnitude ranges, a more explicit expression of volume of damaged rocks (in km3 per yr per km2 of the Earth’s surface at the seismically active region) is and

Evaluating (S10) and (S11) gives V =1.2 10−5 km3 per year per km2 of the Earth’s surface in a seismically active region.

The fraction of lower crust volume affected by earthquake ruptures (1.2%) is obtained using 10 = 4.5, Δε = 5 · 10−5 and the values of other constants mentioned in the text. The result is sensitive to the input parameters, but it is based on realistic values of earthquake intensities in active subduction zones (4.5 annual aftershocks with M > 0 per km2), assumed lower and upper limits of aftershock magnitudes (0 and 8.3), average strain drop (Δε = 5 · 10−5), ratio of damage zone thickness to rupture radius γ = 0.002, and the fraction of aftershocks with hypocenters in the lower crust (1%). Reducing the lower magnitude limit will increase significantly the rupture surface area but not change much the estimated damage zone volume; decreasing the upper magnitude limit from 8.3 to 7.8 will decrease the estimated volume of damaged rock by 32%. Changing the assumed average strain drop Δε = 5 10-5 to average strain drops of 5 10-4 and 5 10-6 will modify the estimated damage volume by factors of 0.21 and 4.64, respectively. Changing γ or the deformation time scale by a given factor (e.g. 2) will modify the damage volume by the same factor. As a consequence, reasonable variations of these parameters will not change significantly the estimated lower crust volume affected by earthquake ruptures.