Interaction between shallow and subcrustal dislocations on a normal fault

Physics of the Earth and Planetary Interiors 129 (2002) 67–82

Interaction between shallow and subcrustal dislocations on a normal fault Andrea Tallarico a,∗ , Michele Dragoni b , Giammaria Zito a a

Dipartimento di Geologia e Geofisica, Università di Bari, Via Edoardo Orabona 4, 70125 Bari, Italy b Dipartimento di Fisica, Università di Bologna, Viale Carlo Berti Pichat 8, 40127 Bologna, Italy Received 8 March 2000; received in revised form 14 May 2001; accepted 17 May 2001

Abstract We propose a model which may explain seismic sequences which are often observed in seismogenic regions, as for example in the Apenninic chain (Italy). In particular, we consider a normal fault and earthquakes taking place at different depths: a first shock in a shallower layer and a second in a deeper one. The normal fault is embedded in a viscoelastic half-space. As a consequence of the rheology, there are two different brittle layers, a shallower and a deeper one, where earthquakes can nucleate. Between these two layers, the rheological behavior is ductile. The thicknesses of the layers depend on the geothermal profile that is calculated taking into account the profile of the thermal and rheological parameters with depth. The fault plane, crossing layers with different rheological behavior, is heterogeneous in respect to the slip style: seismic in the brittle layers, aseismic in the ductile layer. Dislocations in the shallower layer are assumed to produce aseismic slip in the area of the fault belonging to the ductile layer. The stress concentrated, by the seismic and aseismic dislocations, on the fault plane section in the deeper brittle layer is evaluated. It is compared with the tectonic stress rate in order to calculate how much earlier the second earthquake would occur compared to if just the bare tectonic sstress was acting. It results that such an advance is comparable with typical recurrence times of earthquakes and so a mechanism of interaction between different seismic sources, mediated by aseismic slip, can be supposed. The strains and displacements at the Earth’s surface due to seismic and aseismic slip are calculated. They are large enough to be detected by present geodetic techniques. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Seismic sequences; Rheology; Fault dislocation; Stress field; Geotherm

1. Introduction Seismic activity takes place in different ways. Usually seismic shocks are not isolated but happen in sequence and different cases can be recognized. The more frequent case is a main shock followed by aftershocks with very close epicenters concentrated around the asperity that produced the main shock. Aftershocks usually last for periods that are longer, big∗ Corresponding author. Tel.: +39-80-5442627; fax: +39-51-2095058. E-mail address: [email protected] (A. Tallarico).

ger the magnitude of the main shock. Several authors have proposed mechanisms to explain aftershocks and seismic sequences (Richter, 1958). Benioff (1951) showed that in order to explain the temporal aftershock sequence, it was sufficient to assume that rocks suffer a creep-like phenomena and to distinguish between compressional and shearing creep phases. Nur and Booker (1972) proposed an aftershock mechanism based on a time-dependent rock strength, determined by pore pressure diffusing from squeezed pores on the compressive side toward the tensile side of an edge dislocation, thus ascribing the time delay of aftershocks to diffusive processes. Weertman

0031-9201/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 3 1 - 9 2 0 1 ( 0 1 ) 0 0 2 1 1 - 4

68

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

Fig. 1. A sketch of the model, the primed coordinate system lays on the fault plane.

(1974) further refined this model, investigating the interaction of the stress field with an assumed isotropic pore distribution, recognizing the anisotropic influence of crustal permeability on fluid diffusion paths. Bonafede et al. (1980) proposed an aftershock mechanism based on the dilatancy-fluid diffusion theory of the seismic mechanism. Das and Scholz (1981), using the concepts of fracture mechanics, developed a theory of earthquakes which includes delayed multiple events, foreshocks and aftershocks. Hsu et al. (1984–1985) used a two-asperity model to explain the occurrence of aftershocks, supposing migration of a slip front. A mechanism to explain seismic sequences on a strike slip fault was proposed by Dragoni and Tallarico (1992). They considered a mechanism of interaction between seismic sources, placed at the same depth, taking into account the occurrence of aseismic slip. A similar model was proposed by Nomanbhoy and Ruff (1996). The aim of the present paper is to develop a model considering the interaction between different seismic sources placed at different depths on a normal fault in presence of a depth-dependent rheology. A source is placed in a shallower seismogenic layer, the other in a deeper, subcrustal seismogenic layer (Fig. 1). In order to know the position of the

seismogenic layers, we consider the plate boundary embedding the normal fault as a viscoelastic medium able to relax deviatoric stresses (Dragoni et al., 1996; Dragoni et al., 1997; Dragoni et al., 1993). The geothermal profile affects the rheological model through the viscosity. Once evaluated the thickness and depth of the seismogenic layers, we estimate the effect of a shallow dislocation on the deeper seismogenic layer in terms of stress transfer mediated by aseismic slip in the ductile layer. The model considers the surface displacements and strains due to shear rectangular dislocations in a homogeneous half-space.

2. The model In order to apply the rheological model, it is necessary to define the lithospheric geotherm. We assume that the lithosphere is composed of four different layers, having different thermal and/or rheological properties according to their lithological composition. We call zi (i = 1, 2, 3, 4) the depth of the interface between layer i and layer i + 1, while z3 = zc is the base of the crust. The base z4 = z of the lithosphere is defined as the isothermal plane T = Ta = 1330 ◦ C,

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

which is assumed to be the temperature of the asthenosphere. We consider time intervals for which it can be assumed that the geothermal profile T (z) is unchanged. Assuming that in each layer i, the thermal conductivity Ki is uniform, we solve the steady-state heat conduction equation:  H1   − , 0 ≥ z ≥ z1   2 K1 d T = (1)  H0 dz2  −(z1 −z)/D  − e , z1 ≥ z ≥ z K2 with appropriate boundary conditions:  T (0) = 0        T (z+ ) = T (z− )  1 1  dT + dT −   K1 (z1 ) = K2 (z )   dz dz 1      T (z ) = Ta

(2)

where H1 is the radiogenic heat productivity in the layer 1, H0 the radiogenic heat productivity at the top of layer 2 (z = z1 ), and D a length scale for the decrease of radiogenic productivity with depth in layers 2, 3, and 4. We obtain: layer 1 (sedimentary cover with constant radiogenic heat production): H1 2 T (z) = − z − C1 z, 0 ≥ z ≥ z1 (3) 2K1 layers 2, 3, and 4 (remaining upper crust, lower crust, and lithospheric mantle, with radiogenic heat production): H0 2 −(z1 −z)/D D e + C2 (z1 − z) + C3 , K2 (4) z1 ≥ z ≥ z

T (z) = −

C3 =

H0 D 2 −(z1 −z )/D e − (z1 − z )C2 + Ta K2

69

(7)

In order to define the thicknesses and positions of the seismogenic layers, we assume that the lithosphere is subject to a uniform and constant strain rate R and behaves as a viscoelastic medium able to relax the shear stresses (Fig. 2). The exact form of the constitutive equations for the lithosphere is still an open problem. We consider the boundary zone as made of a thermorheologically simple material (Christensen, 1971), i.e. a linear (Maxwell) viscoelastic medium, with viscosity depending on position. This rheology is able to satisfactorily reproduce the gross behavior of the shear zone. A constitutive equation for such a medium is given by Dragoni et al. (1993):   µ 1 t t t t σ˙ ij + δij + 2µe˙ijt (8) σij − σkk δij = λe˙kk η 3 where η is the viscosity, λ and µ the Lamé constants, σijt and eijt , respectively, the stress and strain tensors due tectonic activity, δij the Kronecker delta. Viscosity can be expressed by the relation (e.g. Ranalli, 1987): η(z) = AD exp

E + pV RT

(9)

where E is the activation energy, AD the Dorn parameter, p the pressure, V the activation volume, R the gas constant, and T the absolute temperature. In the lithosphere, one has pV E, so that the term containing the pressure is normally neglected. The solution of the constitutive equation for the extensional case is t t σyy (z, t) = σ¯ yy [1 − 41 (e−t/τ1 + 3e−t/τ2 )]

(10)

where C1 = − C2 =

H1 H0 K2 z1 + D+ C2 K1 K1 K1

K1 z 1 K 2 − z 1 K1 + z  K 1  H1 2 H0 Dz1 × −Ta + z − 2K1 1 K1  H0 D 2 −(z1 −z )/D + (1 − e ) K2

(5)

(6)

Fig. 2. The lithosphere model for an extensional tectonic margin. The dashed line represents the base of the crust.

70

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

where t σ¯ yy (z)

= 4Rη

(11)

t is the asymptotic values of tectonic stress for t → σ¯ yy ∞. The relaxation times appearing in Eq. (10) are η τ1 = (12) µ

τ2 =

λ + 2µ τ1 λ + (2/3)µ

(13)

If we assume λ = µ, we have τ2 = (9/5)τ1 . As a consequence of the discontinuity of E at the base of the crust, the viscosity has a discontinuity at the same depth: so it is smaller in the lower crust than in the lithospheric mantle. Consequently also, the normal t shows a discontinuity at the same depth. stress σyy In order to determine the maximum depth of earthquake nucleation, we employ the differential stress σ (z), difference between the maximum and the minimum principal stresses. We compare the maximum value σ¯ of σ , which is reached for t → ∞, with the minimum value which is necessary to produce slip on a fault with the most favorable orientation (Sibson, 1974): t σ¯ = σ¯ yy

(14)

The value of σ necessary to produce the slip of a generic fault can be expressed as σ ∗ (z) = γ κ(p − p )

(15)

where κ is the friction coefficient, p the pore pressure, and γ a function of the friction coefficient and the fault orientation, while z p(z) = −g ρ(z ) dz (16) 0

γ = ( 21 sin 2δ + κ sin2 δ)−1

(20)

The tangential component of traction on the fault plane is Tfd = Tid mi

(21)

where Tid = σijd nj

(22)

is the traction vector on the fault plane due to the dislocations and n = (0, − sin δ, cos δ)

We assume hydrostatic pore pressure:

where ρ is the density of water.

d σijd = λekk δij + 2µeijd

(23)

(17)

where δ is the dip angle of the fault. The minimum value of σ ∗ is obtained when π 1 1 δ = − arctan (18) 2 2 κ p (z) = −ρ gz

We assume that the coseismic stress due to seismic slip in the shallower seismogenic layer is able to cause or to accelerate aseismic slip on the downward extension of the fault plane, thus concentrating stress in the deeper seismogenic layer. It is possible to evaluate the total stress field due to the seismic and aseismic slip simply by adding the contributions of two dislocations: one is seismic, characterized by a fixed slip U s , and another is aseismic, characterized by a slowly increasing slip U a (t). We consider an elastic homogeneous half-space and calculate, following Okada (1992), displacements udi and strains eijd due to a finite rectangular dislocation on a normal fault. In this way, the effect of the dislocation strain field is evaluated assuming an elastic homogeneous half-space, while the rheological behavior assumed to define the depth extent of the brittle and ductile regions involves a viscoelastic layered half-space. This outward contradiction can be explained taking into account the different scales of time characterizing a seismic sequence and the viscoelastic process given by Eqs. (12) and (13). The stress tensor is easily evaluated using Hooke’s law:

(19)

is the vector perpendicular to the fault plane and m = (0, − cos δ, − sin δ)

(24)

is the vector parallel to the fault plane. Following Eqs. (10), (21) and (22), we can calculate the tangential component of the traction on the fault plane due to the tectonic activity: Tft = Tyt my

(25)

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

Fig. 3. The parameters depending on the lithospheric layers: (a) activation energy E; (b) density ρ; (c) thermal conductivity K.

71

72

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

where

3. Discussion and conclusions

t ny Tyt = σyy

(26)

is the only nonvanishing traction component due to the tectonic activity. The tangential traction on the fault plane Tfd due to coseismic and aseismic slip, calculated by Eq. (21), has to be added to the tectonic traction Tft calculated by Eq. (25). In this way, in the deeper seismogenic layer, there will be a stress concentration increasing with the amount of aseismic slip. For this reason, the effect of the seismic and aseismic slip is to advance the time of a possible second earthquake nucleated in the deeper (subcrustal) seismogenic layer. The advance in time A can be estimated by dividing the traction on the fault plane Tfd , due to both seismic and aseismic slips, by the rate of traction increase on the fault plane T˙ft , due to tectonic stress A=

Tfd T˙ t

(27)

f

where T˙ft (z, t) =

t σ¯ yy

   1 1 −t/τ1 1 −t/τ2 sin δ cos δ e + e 4 τ1 3τ2 (28)

is calculated by Eq. (25).

The present model can be applied to a seismogenic zone characterized by earthquakes taking place on normal faults. A similar condition is present for instance in the Apenninic chain (Italy), where seismicity often occurs in the form of seismic sequences (Amato et al., 1998; Selvaggi and Amato, 1992). The 1997–1998 Umbria-Marche (Italy) earthquake sequence represents a possible application of the present model. It consisted of nine main shocks, the hypocenters of eight were located in the crust at depths of about 5 km, one was located in a subcrustal layer 47 km deep (Cocco et al., 2000); on 9 February 2000, the Italian Seismic Network recorded in the same area, another shock with hypocenter about 50 km deep. Stress transfer between fault segments in the southern Apennines was considered by Nostro et al. (1997) for an elastic medium. In the Apennines extensional stresses act in a zone where, at least in the past, a compressive stress field connected with the uplift of the Apenninic chain was present. In the present model, we use the layer structure shown in Fig. 3, where the activation energy E, the density ρ and the thermal conductivity K are plotted as functions of depth z. The geotherm is shown in Fig. 4 and has been calculated by using Eqs. (3) and (4). The geotherm affects the rheological model through the viscosity η given by Eq. (9), the model sensitivity to different temperature profiles has been investigated by Dragoni (1990). Fig. 5 shows the viscosity as a function of depth. A rheological

Fig. 4. The geotherm. The value of the model parameters D, H1 , H0 , Ta are fixed and shown in Table 1, and the thermal conductivity Ki is shown in Fig. 3c.

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

73

Fig. 5. The viscosity as a function of depth. The value of the model parameter AD is fixed and shown in Table 1, and the activation energy E is shown in Fig. 3a.

t as a function of depth for different times. The dotted line represents the asymptotic stress σ t . The value of the Fig. 6. Normal stress σyy ¯ yy model parameters R, λ, µ are fixed and shown in Table 1.

t (solid line) and the differential stress σ ∗ (dotted Fig. 7. The seismogenic layers obtained by comparison between the normal stress σyy line). The value of the model parameters R, κ, g are fixed and shown in Table 1, and the density ρ is shown in Fig. 3b.

74

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

d; Fig. 8. Components of the stress tensor σijd due to the seismic and aseismic dislocations contributing to traction on the fault plane: (a) σyy d d a (b) σzz ; (c) σzy . Three different values of the aseismic slip U are considered. The black lines represent the intersection of the dislocations with the plane zy. The value of the model parameters W a , W s , La , Ls , λ, µ are fixed and shown in Table 1.

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

75

Fig. 8. (Continued).

Fig. 9. The tangential traction on the fault plane Tfd due to the seismic and aseismic dislocations. Three different values of the aseismic slip U a are considered. The value of the model parameters W a , W s , La , Ls , λ, µ, δ are fixed and shown in Table 1.

76

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

Fig. 10. Rate of tangential traction increase on the fault plane T˙ft as a function of the depth. The value of the model parameters R, λ, µ, δ are fixed and shown in Table 1.

discontinuity appears at the base of the crust: it is due to the discontinuity in the activation energy E. The sharp increase in viscosity at the base of the crust t , given causes a sharp increase of tectonic stress σyy by Eq. (10), as can be seen in Fig. 6 . Fig. 7 shows t and a comparison between the asymptotic stress σ¯ yy ∗ the differential stress σ (z) given by Eqs. (11) and

(15), respectively. It appears that the tectonic stress is greater than the differential stress in two different layers, said brittle, where a dislocation can nucleate. The first layer extends from the free surface to about 15 km, the second from 35 km (the base of the crust) to about 60 km. In the intermediate layer, said ductile, the tectonic stress does not reach the differential

Fig. 11. Advance in time A of the second earthquake on the fault plane x y .

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

77

d ; (b) Fig. 12. The nonvanishing components of the stress tensor eijd due to the seismic and aseismic dislocations on the plane xy: (a) exx d d d a eyy ; (c) ezz ; (d) exy . Three different values of the aseismic slip U are considered. The yellow rectangles and lines represent, respectively, the projections of the dislocations on the plane xy and the fault trace. The value of the model parameters W a , W s , La , Ls , λ, µ are fixed and shown in Table 1.

78

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

Fig. 12. (Continued).

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

79

Fig. 13. The vertical displacement udz at the Earth surface due to the seismic and aseismic dislocations. Three different values of the aseismic slip U a (t) are considered. The yellow rectangles and lines represent, respectively, the projections of the dislocations on the plane xy and the fault trace. The value of the model parameters W a , W s , La , Ls , λ, µ are fixed and shown in Table 1.

80

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

Fig. 13. (Continued).

Fig. 14. The horizontal displacement vector udx + udy due to the seismic and aseismic dislocations. Three different values of the aseismic slip U a are considered. The yellow rectangles and lines represent respectively the projections of the dislocations on the plane xy and the fault trace. The value of the model parameters W a , W s , La , Ls , λ, µ are fixed and shown in Table 1.

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

stress: as a consequence of the viscoelastic behavior, the tectonic stress is relaxed faster than it is accumulated. For this reason, dislocations cannot nucleate in this layer as well as at depths greater than that of the base of the second seismogenic layer. In order to evaluate the contribution to the total stress field due to the dislocations, both seismic and aseismic, we developed a FORTRAN code, implementing the formulas given by Okada (1992), obtaining the strain eijd due to a rectangular dislocation embedded in a homogeneous elastic half-space. The stress σijd is calculated by Eq. (20). For t = 0, we consider a rectangular dislocation, with sides W s = 10 km, Ls = 5 km, dip angle δ = 40◦ , slip U s = 0.5 m placed in the upper seismogenic layer. For t > 0, we add the contribution of a rectangular aseismic dislocation with sides W a = 10 km and La = 31.8 km, taking place on the same fault and extending from the top to the bottom of the ductile layer (Fig. 1). Three different slip amplitudes have been considered. The nonvanishing components of the stress tensor are shown in Fig. 8, where the stress concentration at the top of the lower (subcrustal) seismogenic layer can be seen. The components of the stress tensor can be combined following Eq. (22) obtaining the traction vector on the fault plane. Traction is shown in Fig. 9, where the stress evolution on the fault plane is evident. At t = 0, the seismic dislocation takes place, while the aseismic dislocation is absent and U a (t = 0) = 0, so the first panel of the figure represents traction on the fault plane due to the presence of the seismic dislocation only. The following panels describe the traction on the fault plane for increasing amounts of slip, and in the last panel, the slip of the aseismic dislocation reaches that of the seismic dislocation and the stress transferred to the lower seismogenic layer reaches its maximum value. In this way, a possible earthquake, nucleating in the subcrustal seismogenic layer, will take place earlier compared to if it occurred as a consequence of the bare tectonic stress. The rate of traction increase T˙ft due to tectonic stress, given by Eq. (28), is shown in Fig. 10, where it can be seen that it can be regarded as a constant for depths extending from the Moho to about 70 km. The advance in time A of the lower earthquake, given by Eq. (27), is shown in Fig. 11. It appears that it reaches values of hundreds of years in the subcrustal seismogenic layer; this advance is comparable with typical recurrence times of earthquakes (from tens to hundreds of years). A mech-

81

Table 1 Numerical values of model parameters which are kept constant in drawing the graphs Parameter

Numerical value

Ws

10 10 5 31.8 3×1010 3×1010 9.8 3.15 0.5 1330 8 10−15 2×1011 0.75 40

(km) W a (km) Ls (km) La (km) λ (Pa) µ (Pa) g (m s−2 ) H0 (␮W m−3 ) H1 (␮W m−3 ) Ta (◦ C) D (km) R (s−1 ) AD (Pa s) κ δ (◦ )

anism of interaction between seismic sources placed in the two seismogenic layers, mediated by aseismic slip in the ductile layer, can be supposed. For times comparable with those of a seismic sequence (∼1 year), the behavior is elastic, so the coseismic slip should not experience any important relaxation. This can be inferred by comparison with the characteristic times of the viscoelastic relaxation given by Eqs. (12) and (13) which are in the order of 300 years. In order to check the validity of the model by a comparison with observations, it is useful to evaluate the strain and displacement at the Earth’s surface due to the seismic and aseismic slips. The only nonvanishing components of strain ensor at z = 0, exx , eyy , ezz , and exy , are shown in Fig. 12. Fig. 13 shows the vertical displacement uzd ; in Fig. 14, the horizontal displacement field udx + udy is plotted. With the values employed for various quantities (Table 1), the displacements appear to be of the order of centimeters, and so they can be detected by current geodetic techniques.

Acknowledgements This research was funded by the Italian Space Agency under grant ARS-98-68 and by MURST. Many thanks are due to Francesco Mongelli for useful comments and suggestions on the first version of the manuscript.

82

A. Tallarico et al. / Physics of the Earth and Planetary Interiors 129 (2002) 67–82

References Amato, A., Azzara, R., Chiarabba, C., Cimini, G.B., Cocco, M., Di Bona, M., Margheriti, L., Mazza, S., Mele, F., Selvaggi, G., Basili, A., Boschi, E., Corboulex, F., Deshamps, A., Gaffet, S., Bittarelli, G., Chiaraluce, L., Piccinini, D., Ripepe, M., 1998. The 1997 Umbria-Marche, Italy, earthquake sequence: a first look at the main shocks and aftershocks. Geophys. Res. Lett. 25, 2861–2864. Benioff, H., 1951. Earthquakes and rock creep, Part 1. Bull. Seism. Soc. Am. 41, 31–62. Bonafede, M., Mulargia, F., Boschi, E., 1980. Implications of the dilatancy-fluid diffusion theory for aftershocks sequences. Il Nuovo Cimento 3, 180–190. Christensen, R.M., 1971. Theory of Viscoelasticity: An Introduction. Academic Press, New York. Cocco, M., Nostro, C., Ekström, G., 2000. Static stress changes and fault interaction during the 1997 Umbria-Marche earthquake sequence. J. Seismol. 4, 501–516. Das, S., Scholz, C.H., 1981. Theory of time-dependent rupture in the Earth. J. Geophys. Res. 86, 6039–6051. Dragoni, M., 1990. Stress relaxation at the lower dislocation edge of great shallow earthquakes. Tectonophysics 179, 113– 119. Dragoni, M., Tallarico, A., 1992. Interaction between seismic and aseismic slips along a transcurrent fault: a model for seismic sequences. Phys. Earth Planet. Int. 72, 49–57. Dragoni, M., Santini, S., Tallarico, A., 1993. A viscoelastic shear zone model of compressional and extensional plate boundaries. Pure Appl. Geophys. 140, 471–491.

Dragoni, M., Doglioni, C., Mongelli, F., Zito, G., 1996. Evaluation of stresses in two geodynamically different areas: stable foreland and extensional back-arc. Pure Appl. Geophys. 146, 319–341. Dragoni, M., Harabaglia, P., Mongelli, F., 1997. Stress field at a transcurrent plate boundary in the presence of frictional heat production at depth. Pure Appl. Geophys. 150, 181–201. Hsu, V., Helsley, C.E., Berg, E., Novelo-Casanova, D.A., 1984– 1985. Correlation of foreshocks and aftershocks and asperities. Pure Appl. Geophys. 122, 878–893. Nomanbhoy, N., Ruff, L.J., 1996. A simple discrete element model for large multiplet earthquakes. J. Geophys. Res. 101, 5707– 5723. Nostro, C., Cocco, M., Belardinelli, M.E., 1997. Static stress changes in extensional regimes: an application to southern Apennines (Italy). Bull. Seism. Soc. Am. 87, 234–248. Nur, A., Booker, J.R., 1972. Aftershocks caused by pore fluid pressure. Science 175, 885–887. Okada, Y., 1992. Internal deformation due to shear and tensile faults in a half-space. Bull. Seism. Soc. Am. 82, 1018–1040. Ranalli, G., 1987. Rheology of the Earth. Allen & Unwin, Boston. Richter, C., 1958. Elementary Seismology. Freeman, San Francisco, CA. Selvaggi, G., Amato, A., 1992. Subcrustal earthquakes in the northern Apennines (Italy): evidence for a still active subduction. Geophys. Res. Lett. 19, 2127–2130. Sibson, R.H., 1974. Frictional constraints on thrust, wrench and normal faults. Nature 249, 542–544. Weertman, J., 1974. Water flow around an earthquake fault dislocation. J. Geophys. Res. 79, 2391–2393.