Strain rate relaxation of normal and thrust faults in Italy

S U M M A R Y We find that geodetic strain rate (SR) integrated with the knowledge of active faults points out that hazardous seismic areas are those with lower SR, where active faults are possibly approaching the end of seismic cycle. SR values estimated from GPS velocities at epicentral areas of large historical earthquakes in Italy decrease with increasing elapsed time, thus highlighting faults more prone to reactivation. We have modelled an exponential decrease relationship between SR and the time elapsed since the last largest earthquake, differencing historical earthquakes according to their fault rupture style. Then, we have estimated the characteristic times of relaxation by a non-linear inversion, showing that events with thrust mechanism exhibit a characteristic time (∼ 990 yr) about three times larger than those with normal mechanism. Assuming standard rigidity and viscosity values we can infer an average recurrence time of about 600 yr for normal faults and about 2000 yr for thrust faults.


I N T RO D U C T I O N
The recent significant earthquakes occurred in Italy at L'Aquila (2009 April 6, 13.38E 42.34N, Mw = 6.3, central Apennines, normal mechanism) and in Emilia region (2012 May 20, 11.23E 44.89N, Mw = 6.1, northern Apennines, thrust mechanism) are the striking consequence of Adriatic subduction and rollback (e.g. Devoti et al. 2008;Carminati et al. 2012). The slab retreat generates compression along the frontal Apennines accretionary prism, moving towards NE at a rate of 2-3 mm yr -1 Bennett et al. 2012), and extension along the Apennines chain (Doglioni 1991). All these mechanisms generate an inhomogeneous deformation belt, characterized by some relative maxima and minima in the strain rate (SR) field, regardless the extensional or compressive tectonic style (Fig. 1). The along-strike variations in the accretionary prism are controlled by the heterogeneities of the inherited Tethyan passive continental margin, determining short-wavelength tectonic features, hence preventing large magnitude events.
Both the above-mentioned earthquakes occurred in areas of lower SR with respect to the neighbourhood. In fact, GPS networks allow to estimate the SR of ∼19·10 −9 yr −1 around L'Aquila and ∼14·10 −9 yr −1 in Emilia, in the epicentral area during the preseismic period, both almost above the significance level.
Several geological and geodetic observations have been made to model the seismic cycle, recording the coseismic and interseismic crustal deformation consequent to rupture and locking stage of faults. GPS measures the displacements within the earthquake cycle, but due to the viscoelastic behaviour of the lower crust and mantle lithosphere, the surface deformations vary throughout the seismic cycle. Early in the cycle, viscous flow in the lower crust exerts a traction on the upper crust and produces high-velocity gradients across active faults, that is, high SR. Late in the cycle the opposite occurs and SR decreases across the faults (Nur & Mavko 1974;Savage & Prescott 1978). Remarkable achievements have been accomplished in studying time-varying SRs over transform faults (e.g. Thatcher 1983;Scholz 1992;Dixon et al. 2002;Hetland & Hager 2006;Hilley et al. 2009;Lundgren et al. 2009;Barbot et al. 2012), but only a limited number of papers deals specifically with changes of the deformation field induced by thrust faults (Thatcher & Rundle 1979;Wang et al. 2003;Pollitz et al. 2010) and none, in our knowledge, concerning normal fault systems. Faure Walker et al. (2010) found significant discrepancies between geodetic and geological observed SRs in the central Apennines and coherently called for a temporal variability in the rate of interseismic strain accumulation.
Purely elastic models used to study interseismic deformation and to analyse hazard potential for future earthquakes are mainly based on the implicit assumption that interseismic deformation rate does not change with time (Wang 2007). Some key observations conducted in highly active subduction areas, for example, in the Cascadia region (Dragert et al. 1994) and in Peruvian Andes (Perfettini et al. 2005), give evidence that the deformation rate is strongly time-dependent, so that geodetic observations made at a given time, represent a 'snapshot' of an evolving deformation field. Hammond et al. (2010) determined the presence of a long-term relaxation signal superimposed on the background secular velocity field Figure 1. SR map of the Italian area (contour lines every 20·10 −9 yr −1 ) obtained from the 2D velocity field . The grey line with black triangles marks the active fronts of deformation, delimiting the 'undeformed' foreland of the Adriatic microplate in the footwall. The Adriatic microplate subducts both W-wards, beneath the Apennines, and E-wards under the Dinarides; it overrides NW-wards the Eurasian Plate to form the Alps. Yellow stars indicate the locations of the recent 2009 April 6 L'Aquila and 2012 May 20 Emilia events. The black boxes are the individual seismogenetic sources reported in the DISS database (Basili et al. 2008). The last significant seismic events reported in CPT04 (Rovida et al. 2011) generated by the selected faults of the DISS database are the blue (normal faults) and the red (thrusts) boxes. estimated from GPS time-series. Wang et al. (2001) demonstrated that geodetic observations in the interseismic period can be modelled introducing viscoelasticity or an exponential decreasing correction in the downdip transition zone (Wang et al. 2003), since purely elastic dislocation model cannot account for all the slip deficit loaded by plate convergence rate. The most important effect of introducing viscoelasticity in the model is the time dependence of the interseismic deformation, so that the SR is expected to decrease with time (Wang et al. 2001). At fault level, the brittle-ductile transition plays as a switch in controlling the seismic cycle; the steady-state shearing in the ductile lower crust determines the loading of elastic energy in the volume near the locked section of the fault, in the brittle upper crust, where the SR is relatively lower (Doglioni et al. 2011).
Since a large earthquake may rupture only a segment of a large tectonic structure, different segments along the same margin may be at different stages of strain accumulation towards future earthquakes, as shown in the case of the interseismic loading area detected before the Chilean Maule earthquake (Madariaga et al. 2010).
Therefore, the SR mapping from interseismic GPS velocities in tectonically active areas may be considered a possible indicator of seismic gaps, that is, areas where the tectonic loading is active and faults are presumably locked and loading elastic energy late in the seismic cycle. Such idea has been conveniently applied in the Italian area to point out areas of possible future large seismic events (Riguzzi et al. 2012).
Since the interseismic SR values based on the last decade of GPS observations, are 'snapshots' of a long term, possibly nonlinear, trend throughout the seismic cycle, they could be used to infer the time evolution of deformation in different segments of the seismogenic structures. Anyway, since SR lows can be found also in regions where no present tectonics is expected, the GPS-derived by guest on August 23, 2013 http://gji.oxfordjournals.org/ Downloaded from SR field must be necessarily integrated with the present knowledge of active faults (regardless they reach the surface or are buried blind faults).
Geodetic measurements devoted to active tectonics studies are extensively carried out in Italy only since the early 2000 s, the maximum effort in deploying the GPS networks in Italy dating back to about 2005 (Avallone et al. 2010), so that the observed deformations represent the instantaneous seismic cycle conditions in every area of the Italian region (e.g. Devoti et al. 2012).

M E T H O D
To characterize the long-time behaviour of interseismic SR, we computed the SR 'snapshot' over selected faults that correspond to the epicentral areas of large historical earthquakes with M ≥ 5.6. We have chosen only those seismogenic structures (Basili et al. 2008), where no seismic event of 5.0 ≤ M < 5.6 occurred since the epoch of the large main shock (M ≥ 5.6). In our experience, M = 5.0 can be considered as the lower limit over which coseismic deformations could be recorded by GPS, and thus able to perturb SR.
At each 'instantaneous' value of SR we have associated the elapsed time of the last significant earthquake (as in Thatcher 1993), separating the normal faults from the thrust planes, to study the obtained time trend. The underlying hypothesis is that the interseismic SR on different seismogenic sources of the same type behave as a stationary process (i.e. their time decaying rate are similar), dependent only on the time elapsed since the last strong earthquake. Therefore, we can estimate the time decaying model parameters just sampling the process on different seismogenic sources for which different intervals since the last strong earthquake elapsed.
Such processes can be described by a non-linear decay law (Savage & Prescott 1978), corresponding to constitutive equations of viscoelastic materials, where we have tested a classical process of time exponential decrease, separating the normal from thrust faults.
We have described the normal and thrust fault time evolution with an exponential decay law S R(t) = a · e −b·t , typical of processes in which SR decreases at a rate proportional to SR itself, where τ = b −1 can be considered as characteristic time of the process for which the SR assumes the value S R(t = τ ) = a e , and a the asymptotic value of SR at infinite. The model is characterized by two parameters a = (a n , a t ) and b = (b n , b t ) according to the fault style (n = normal or t = thrusts).

DATA A N D M O D E L
We have computed the 2-D strain rate tensor of the Italian area from the planar GPS velocities and the associated uncertainties ) by a distance-weighted approach using all stations on a regularly spaced grid downweighting with the function W = exp(−d 2 /α 2 ) that depends on the spatial distribution of GPS sites (Shen et al. 2007). The average inland smoothing distance applied to the Italian Peninsula was about 50 km. We have computed the deformation rate (SR) as the scalar S R = (ε 2 11 +ε 2 22 + 2 ·ε 2 12 ) accounting for all the 2-D deformation tensor components, thus representing the total amount of deformation rate. In the Italian area, SR ranges from zero to about 90 · 10 −9 yr −1 with maximum average values along the Apennine belt, around the Messina Strait, the NE Alpine sector and Dinarides (Fig. 1). The SR average uncertainty is about 9 · 10 −9 yr −1 .
Therefore, we have interpolated the value of SR at the epicentral coordinates of some large historical earthquakes. The Italian area is characterized by moderate seismicity with moment magnitudes never exceeding 7.4. Thanks to the good record of historical and instrumental seismicity, we have selected from the CPT04 catalogue (e.g. Rovida et al. 2011) all the seismic events matching the following criteria: (1) occurred after 1000 a. C. with magnitude Mw ≥ 5.6 (191 earthquakes) on well-recognized faults reported in the DISS database (Basili et al. 2008); (2) the faults should have a pure normal or thrust (reverse) tectonic style; (3) each event should never be followed by significant events (Mw ≥ 5.0), so to be considered the last significant fault event; (4) offshore and deep earthquakes, events occurred in volcanic areas and near the border of the SR map of Fig. 1 were discarded.
Each event is then assigned to a corresponding fault characterized by a pure extensional or contractional mechanism activated in the past, not experiencing recent significant seismic activity.
The seismic events survived after our severe selection are reported in the Table 1. We added to the list, the 1570 Ferrara event (Mw = 5.5) occurred near the epicentral area of the recent Emilia seismic sequence (2012 May), since of recent interest. Fig. 1 shows the SR map of the Italian area, the black boxes are the individual seismogenetic sources of the Database of Individual Seismogenic Sources (DISS); the selected seismic events are the red (thrusts) and blue (normal faults) boxes.
The time elapsed from each event to now is then associated to SR at its epicentral area: both the normal (blue squares) and thrust mechanism (red squares) events of Fig. 1 clearly show a decreasing SR with increasing elapsed time, with different decaying time. Fig. 2 shows the SR sampled at the epicentres of the seismic events (see the Table 1) versus the elapsed time (blue bullet, normal faults and red diamond, thrusts).
) 2 as indicator of the goodness of fits, where N is the number of data, p the number of estimated parameters, S R mod i and S R obs i the modelled and observed SR values, σ i the observation errors. The χ 2 values of the two least-squares fits are 3 and 0.8, respectively, for normal and thrust case. The SR values at t = τ are S R n (τ n ) = 20.1 · 10 −9 yr −1 and S R t (τ t ) = 8.5 · 10 −9 yr −1 , the first matches the more frequent value of SR in the Italian area (Riguzzi et al. 2012), the second is around the average significance level.

D I S C U S S I O N
According to classical viscoelastic models describing a simple earthquake cycle (pre-earthquake strain accumulation, coseismic strain release and post-seismic adjustment), SR is expected to decrease with time during the interseismic period, nevertheless, a complete relaxation is never reached because the process is asymptotic and the recurrence time of earthquakes may be less than the relaxation time of the material (e.g. Savage & Prescott 1978;Wang et al. 2001).  (Table 1) and elapsed time (blue bullet, normal faults and red diamond, thrusts). The model interpolating the observations is an exponential decreasing function with different decay rates. The characteristic times are τ n = 317 ± 204 and τ t = 991 ± 690 yr, respectively, for normal faults and thrusts. Therefore, thrusts have a relaxation time about three times larger than normal faults. Laboratory experiment and natural faults have shown that static friction and healing vary with loading rate and time, and the postseismic healing is expected to be retarded for a long time after an earthquake (Marone 1998). Recent geodetic observations in the North American continental interior give evidence that fault loading, strength, or both vary with time (Calais & Stein 2009). Caporali et al. (2011) provided indications of a weak linear correlation between the geodetic shear strain rate and the time elapsed from the last largest earthquake; however, their analysis was performed on 36 extensive areas (the seismic zones depicted in Meletti et al. 2008), rather than isolated faults. On the contrary, according to Riguzzi et al. (2012), the spatial pattern of SR estimated from interseismic GPS velocities can be considered a useful indicator of areas where the tectonic loading is active and faults are late in the seismic cycle, towards a locking state. Then, we have estimated the interseismic SR from GPS velocities in the Italian area, considering its spatial pattern just like 'snapshots' of a long-term non-linear time trend during the seismic cycles.
We have modelled an exponential decrease relationship between SR and the time elapsed since the last largest earthquake for normal and thrust faults, estimating the characteristic times of relaxation by a non-linear inversion. The events with thrust mechanism exhibit a characteristic time about three times larger, and lower SR, than those with normal mechanism. Such differentiation is explained by the need to reach a larger amount of tectonic loading to activate ruptures in compressional environments with respect to the energy budget required to generate ruptures by traction. This evidence is congruent with the conclusion inferred in Kanamori & Allen (1986) and Romanowicz (1983) about the involvement of larger energy amounts and longer times to activate thrusts with respect to normal faulting processes.
We hypothesize that our characteristic times are not so far from the repeat time of earthquakes, in fact according to classical by guest on August 23, 2013 http://gji.oxfordjournals.org/ Downloaded from viscoelastic models the repeat time of earthquakes is controlled by the rate of tectonic loading and non-linear stress relaxation mechanisms, related to the characteristic time. Perfettini & Avouac (2004) showed that a periodic seismic cycle depends essentially on the fault frictional properties, afterslip and viscous relaxation and the repeat time is proportional to viscosity. The introduction of the time dependence of deformation following large earthquakes and palaeoseismic estimates of fault slip rates in a Maxwell viscoelastic model allowed Segall (2002) to estimate the recurrence interval of about 280 yr along the San Andreas Fault.
According to a simple model of stick-slip motion of an elastic plate (the lithosphere) overlying a viscoelastic substrate (Savage & Prescott 1978), the more the depth of seismogenic zone is comparable to the thickness of the elastic layer, the more significant is the effect of relaxation. Following this model and associating the characteristic time to the relaxation time of lower crust, it is possible to infer the viscosity of the substrate as η = τ ·μ 2 , where τ is the relaxation time and μ the rigidity. Savage & Prescott (1978) have shown that the recurrence time of earthquakes is proportional to the relaxation time of the lower crust, for example, applying their model to sections of the San Andreas Fault, the ratio is 2.4 if η is 4 · 10 19 Pa · s.
Consequently, extending the reasoning to the Italian case, it is possible to infer that the recurrence time of thrusts is about three times larger than those of normal faults, the same ratio revealed by relaxation times.
Close constraints to the recurrence times would be given by realistic values of rigidity and viscosity. Assuming a standard μ of 3 · 10 10 Pa and an average viscosity value of 10 19 Pa · s (Gardi et al. 2003;Dalla Via et al. 2005;Riva et al. 2007;Ismail-Zadeh et al. 2010), we can infer a recurrence time of about 600 yr for normal faults and approximately 2000 yr for thrust faults.
Of course, our speculations and the idealization implicit in fitting model curves to data in Fig. 2 are worth emphasizing. The modelling results hold only if the data are taken from different periods of identical earthquake cycles characterized by a similar recurrence interval. Because GPS record in Italy is very short (15 yr or less) compared with the high variable earthquake recurrence times inferred from palaeoseismology, geology and statistical seismology (100-10 000 yr) (Pace et al. 2006;Galli et al. 2008;Akinci et al. 2009;Faure Walker et al. 2010, and references therein), the longterm SR profile predicted by classical models is difficult to be constrained. However, if one assumes that the earthquake cycles on the normal and thrust faults are basically similar to a stationary process, then we believe that data from each style can be combined to sample most of the cycle.

C O N C L U S I O N S
We have found and modelled an exponential decrease relationship between SR and the time elapsed since the last largest earthquake for normal and thrust faults in the Italian area. We have estimated the characteristic times of relaxation by a non-linear inversion showing that events with thrust mechanism exhibit a characteristic time about three times larger, and lower SR, than those with normal fault mechanism. We conclude that GPS SR and the knowledge of active faults appear as a useful tool to point out hazardous seismic areas as those with lower SR, where active faults are possibly approaching the end of seismic cycle. Even if any deduction about earthquake repeat time appears at present unripe, we consider it worthwhile to analyse the time-varying SR as an important step towards a better definition of the seismic hazard map of the Italian region.

A C K N O W L E D G E M E N T S
We are grateful to all the people working on the seismological databanks and RING GPS network operated by Istituto Nazionale di Geofisica e Vulcanologia. The editor J. Trampert and an anonymous reviewer are warmly acknowledged.