Viscoelastic relaxation and long-lasting after-slip following the 1997 Umbria-Marche (Central Italy) earthquakes

Summary

We combine Global Positioning System (GPS) measurements with forward modelling of viscoelastic relaxation and after-slip to study the post-seismic deformation of the 1997 Umbria-Marche (Central Apennines) moderate shallow earthquake sequence. Campaign GPS measurements spanning the time period 1999–2003 are depicting a clear post-seismic deformation signal. Our results favour a normal faulting rupture model where most of the slip is located in the lower part of the seismogenic upper crust, consistent with the rupture models obtained from the inversion of strong motion data. The preferred rheological model, obtained from viscoelastic relaxation modelling, consists of an elastic upper crust, underlain by a transition zone with a viscosity of 10¹⁸ Pa s, while the rheology of deeper layers is not relevant for the observed time-span. Shallow fault creep and after-slip at the base of the seismogenic upper crust are the first order processes behind the observed post-seismic deformation. The deep after-slip, below the fault zone at about 8 km depth, acting as a basal shear through localized time-dependent deformation, identifies a rheological discontinuity decoupling the seismogenic upper crust from the low-viscosity transition zone.

crustal deformation, earthquakes, Global Positioning System (GPS), rheology

1 Introduction

The 1997 Umbria-Marche seismic events represent the strongest earthquake sequence in Italy in the last two decades.

The Umbria-Marche region, part of the Central Apennines, is within a complex plate boundary between the Tyrrhenian and the Adriatic in Central Italy. The regional structural setting is characterized by the presence of normal faults cutting the upper crust, related to the backarc extension of the west directed Adriatic subduction with a compressional front progressively migrating eastward (Pialli et al. 1998; Doglioni et al. 1999). The post-orogenic extension has been accompanied by thinning of the crust, volcanic activity and high heat flow (Miocene to Plio-Pleistocene) towards the Tyrrhenian (Pialli et al. 1998, and references therein). In Plio-Pleistocene, extension affected the Apennines, where graben-like structures trending NW–SE overprint the compressional features. On the Adriatic margin, in contrast, the outer front of the belt still exhibits thrust faulting (Meletti et al. 2000).

Two moderate size crustal earthquakes struck the Umbria-Marche area on September 26: the first event (M_w 5.7, 0:33 UTC) was followed by a nearby second large shock (M_w 6.0, 9:40 UTC). A third major event took place on October 14 (M_w 5.6) at 15 km SE of the September events. All the mentioned events, together with several aftershocks (including three with M_w ≥ 5.0), exhibited NW trending normal fault focal mechanisms and relatively shallow hypocentral depths (less than 7 km) (Amato et al. 1998; Chimera et al. 2003). Focal mechanisms and epicentral locations of the three largest events are displayed in Fig. 1.

Figure 1.

GPS network, focal mechanisms for the three main earthquakes (from Capuano et al. 2000) and topography of the area. Rectangles represent the surface projections of the faults, dipping SW, as described in Section 3 (solid blue line for Zollo FM, dashed red line for Salvi FM). Campaign GPS sites are indicated by circles and diamonds, respectively belonging to the first transect (Tr.A) or the second transect (Tr.B). The continuous station of CAME is indicated by a full square.

In 1999 we set-up a first network of six campaign Global Positioning System (GPS) sites distributed along a 30-km-long transect crossing the fault responsible of the largest normal faulting earthquake (1997 September 26, 9:40, M_w 6.0). This network was first occupied in 1999 October and in 2000 September.

The results of the first two campaigns were published by Aoudia et al. (2003). The authors showed that relaxation in the crust explains part of the observed post-seismic deformation. They argued that the localized deformation may require other processes such as after-slip or poro-elastic relaxation.

The existing network was extended in 2000, with the addition of a second transect of four sites crossing the neighbouring fault, responsible of the first event of 1997 September 26 (0:33, M_w 5.7). The network extension was preceded by the installation of a nearby continuous GPS station at Camerino (CAME), as a support to campaign measurements. During the 2000 September campaign, we occupied the network extension for the first time; later, we remeasured the whole network in 2001 September and 2003 May.

In this paper, we extend the results of Aoudia et al. (2003) by means of a larger geodetic data set and a refined modelling of the post-seismic deformation. For the geodetic measurements we extend the GPS data set with the addition of two GPS campaigns and the continuous GPS station of CAME. Therefore, we move a step forward with respect to the 1-D baseline approach adopted by Aoudia et al. (2003) and we describe deformations in terms of 2-D (planimetric) displacements. For the modelling of post-seismic deformation, further to viscoelastic relaxation, we explore the contribution of after-slip and poro-elastic processes.

We test the statistical significance of the measurements and compare them with the predictions of several deformation models. Accordingly, we propose a viscosity structure for the crust and discuss different published coseismic fault rupture models.

2 Seismotectonic Setting

The main active faults of the Umbria-Marche fault zone have been reported by the GNDT working group (Galadini et al. 2000), and those recognized as seismogenic are extending along a NW–SE trend, in an echelon distribution, from Gualdo Tadino to the north up to Norcia to the south, crossing the Colfiorito basin where the large 1997 September 26 events are located. These faults, dipping towards the southwest, define a half-graben structure and form a series of intermontane basins, which alternate with NW–SE trending ranges.

At the southernmost edge of the fault system that has been reactivated during the 1997 earthquake sequence, Meghraoui et al. (1999) reported NE–SW trending fold axes. The Umbria-Marche fold and thrust belt is, therefore, likely to be segmented into three main structural bodies that could explain the interplay between the last three moderate earthquake sequences in the region: Colfiorito in 1997, Norcia to the south in 1979 (MS = 5.8; Deschamps et al. 1984) and Gubbio to the north in 1984 (MS = 5.2; Haessler et al. 1988).

At the scale of the 1997 reactivated fault system, the moment tensor inversions of the two large September 26 and October 14 events show dominantly normal faulting mechanisms, whereas selected aftershocks (magnitude in the range between 2.7 and 4.4) within the Colfiorito basin, reveal that the prevailing deformation at the step-over is of strike-slip faulting type (Chimera et al. 2003).

According to Cello et al. (2000), this step-over zone is marked by pre-existing transverse faults. Furthermore, within the same area, Cinti et al. (2000) report several differently oriented cracks interpreted as the surface effect of minor displacements along transverse structures that are likely to be oriented N–S and may correspond to western edge of the Colfiorito basin. Therefore, it is likely that the step-over between the 1997 reactivated fault fragments (Meghraoui et al. 1999), collocated between the rupture areas of the two September 26 events, is of strike-slip type and could have controlled the lateral propagation of slip.

Shear wave velocity sections across the zone of the 1997 normal faulting earthquake sequence show that the reactivated SW dipping fault zone, delineated by the earthquake foci of the September earthquakes, displays a typical thrust fault geometry (Chimera et al. 2003). Therefore, the 1997 Umbria-Marche normal fault zone can be interpreted as an inversion of pre-existing thrust faults, where structural changes in the upper crust seem to control the present fault characteristics such as: rupture geometry, pattern of deformation and emergence towards the surface.

3 Fault Models

A number of papers have been published based on regional and local seismic networks that estimate the source parameters and discuss the tectonic significance of the 1997–1998 earthquake sequence (e.g. Amato et al. 1998; Ekström et al. 1998; Meghraoui et al. 1999; Barba & Basili 2000; Amato & Cocco 2000, and papers therein; Chiarabba & Amato 2003; Chiaraluce et al. 2003; Chimera et al. 2003).

Other papers provide fault rupture models using strong motion and geodetic data sets. Zollo et al. (1999) propose a rupture model for the two events on 1997 September 29 by inverting strong motion data recorded at near-source distances, later updated by Capuano et al. (2000) with the addition of the 1997 October 14 earthquake. A number of authors (e.g. Hunstad et al. 1999; Stramondo et al. 1999; Salvi et al. 2000; Lundgren & Stramondo 2002; Belardinelli et al. 2003; Crippa et al. 2006) used GPS and DInSAR data. Basili & Meghraoui (2001) and De Martini et al. (2003) use levelling profiles.

The noticeable differences in the various fault rupture models (Table 1) stands between the strong motion and the geodetic models. Therefore, we decided to make use of the results published by Salvi et al. (2000) as far as a geodetic model is concerned (hereafter defined as Salvi FM), while we took the models published by Zollo et al. (1999) and Capuano et al. (2000) as reference seismological model (hereafter defined as Zollo FM). The refinement of Zollo's model proposed by Basili & Meghraoui (2001) (hereafter defined as Zollo-Basili FM), consisting in a high-angle updip extension of the rupture (Table 1), is also considered.

Table 1.

Fault rupture models for the three largest earthquakes.

4 Gps Network And Measurements

Fig. 1 shows the GPS network covering the study area. The monumentation of each site was carefully performed to ensure a submillimetre centring. At each station point, a 25 cm long still rod was fixed in solid rock and a properly designed steel pillar holding the antenna was centred on this ground part. During all campaigns, data were collected over a period of four consecutive days, with daily sessions of 8 hr, contemporaneously for all sites. Data analysis was performed with the Bernese software (Hugenobler et al. 2004), where the Quasi Iono Free strategy for ambiguity fixing was selected, among the other possible options. Due to the high ionospheric activity, mainly in 2001 and 2002, global ionospheric models estimated by CODE (Hugenobler et al. 2000) were used in the L1 and L2 ambiguity estimation step. Tropospheric parameters were estimated on a two-hourly basis and wet zenith delays were modelled as stochastic parameters. Observations were first analysed on a daily basis with multibase approach, yielding full network solutions; subsequently, the daily network normal equations were combined into multiday solutions for each year using the ADDNEQ program of the Bernese software. In order to have realistic errors, formal sigmas, intrinsically small, have been multiplied by a factor ∼3 as inferred from the analysis of the daily coordinate solutions repeatability (Table 2).

Table 2.

Formal and re-estimated errors (mm) of the GPS solutions.

4.1 Antenna phase-centre variations and reference frames

In Aoudia et al. (2003) the results of the first two GPS campaigns for the first transect were presented. The measurements at all sites were in good agreement with model predictions, except for site MONT. The anomalous behaviour of MONT revealed to be caused by antenna mixing during the first two campaigns. Only at MONT, in fact, two different antennas have been used in 1999 and in 2000, namely a TRM22020.00+GP and a TRM33429.00+GP, respectively. All the other sites have always been measured with the same antennas.

In order to verify if antenna mixing could justify the behaviour of MONT, proper tests have been performed by Barzaghi & Borghi (2004). The results of these tests show that antenna mixing using relative Phase-Centre Variation (PCV) parameters is critical in terms of repeatability. By processing the data with relative PCV parameters, differences in the height component up to 7 mm were estimated. On the contrary, submillimetre coordinates differences were obtained by processing the same data with the absolute PCV parameters provided by the German company GEO+ + GmbH (Menge et al. 1998).

The reprocessing of the GPS data published by Aoudia et al. (2003) using these absolute PCV calibration led to an estimate in the MONT baseline variation which amounts to 5 mm (instead of 8 mm obtained with the relative parameters provided by the US National Geodetic Survey). This new result is in better agreement with the model predictions of Aoudia et al. (2003).

Consequently, all the GPS campaigns (1999 October and 2000 September) and the new campaigns (2001 September and 2003 May) have been processed (or reprocessed) using the absolute GEO+ + PCV parameters. Furthermore a new continuous station, established in Camerino (CAME), has been added to the 10 points of our campaign network. CAME has been set up in 2000 April.

All the campaigns have been framed in the ITRF97 and ITRF00 (Altamimi et al. 2002) (ITRF97 for the first three campaigns and ITRF00 for the last one, respectively), using the precise ephemerides provided by the International GNSS Service (IGS) and the IGS continuous stations MEDI and UNPG. Transformation parameters between frames (ftp://lareg.ensg.ign.fr/pub/itrf/ITRF.TP) have been applied to allow the comparison between all campaigns. However, as we are dealing with a local network, any slight mismodelling in the transformation parameters has a negligible impact on relative displacements (computed with respect to CAME).

In fact, if we consider the general transformation formula between frames

1

where are the coordinates into two different frames (in our case, ITRF97 and ITRF00, respectively); is the translation vector computed at the actual epoch, taking into account its rate from the reference epoch 1997.0; δR is the rotation matrix for small rotation angles (r₁, r₂ and r₃) and is given by

2

D is the scaling factor; we have, for relative displacements,

3

As the rotation angles r₁, r₂ and r₃ between the ITRF00 and ITRF97 and their rates are very small (the only non-zero term is the rate of r₃ with a value of 0.02 millisecondarc yr⁻¹), the term is negligible because in our network the maximum distance is about 39 km (baseline SPEL—CAME). The same holds for the term , as the factor (D) and its rate are 1.55 and 0.01 ppb, respectively. So, the relative displacement vector has, at first order, the same components in the two different reference frames.

4.2 Displacement vectors for the whole network

Due to the availability of a larger network and of the CAME continuous station for the campaigns starting from year 2000, we decided to proceed with vector displacements. In this way, even if we neglect the 1999 campaign, we gain a realistic description of the post-seismic deformation process. The baseline approach, in fact, has major shortcomings: it neglects the component of motion perpendicular to the baseline and is intrinsically biased by the low accuracy of the vertical component. Nonetheless, for sake of completeness, at the end of Section 7 we shortly show and discuss rates of post-seismic deformation along a profile.

Campaign results are displayed in Fig. 2, separately for years 2001–2000, 2003–2001 and 2003–2000.

Figure 2.

GPS displacements with respect to CAME. Error ellipses represent 1σ. Surface fault projections are indicated with solid line rectangles for Zollo FM and with dashed line rectangles for Salvi FM.

We see from Fig. 2 how the far-fault sites show a consistent motion through the years, while the near fault sites reveal a different behaviour. The near fault sites CERE, COLL and MONT show an inversion of the direction of motion between 2001–2000 and 2003–2001. The near fault sites on the second transect, POPO, DIGN and CENT, exhibit an almost 90° rotation between 2001–2000 and 2003–2001. RASI is the only site not showing any motion between 2001 and 2003.

The availability of uniformly measured and processed GPS data clearly highlights a long-lasting control of the deformation localized around the fault, superimposed to a longer wavelength deformation process, as exhibited by the difference in motion of near fault sites (CERE, COLL, MONT and POPO, DIGN and CENT) when compared to far-fault sites (SPEL, VALL and SEFR).

4.3 Statistical comparison between geodetic and model deformations

In order to compare the deformations from geodetic measurements and those predicted using different geophysical models, we have applied two different test statistics. The rationale for an accurate characterization of the statistical tests, besides our belief that statistics represents a powerful tool for the interpretation of any set of measurements, lies in the fact that the magnitude of the displacements that we are observing is challenging the capabilities of GPS campaigns. Therefore, we have opted for a very accurate statistical approach to the analysis of our data set, aiming to extract all the information it might contain.

In the first case, after defining the vector as the GPS horizontal displacements in the (N,E) coordinates between two epochs, that is, , and the corresponding coming from a geophysical model, we want to test the hypothesis

4

Assuming that the estimates are normally distributed, it can be shown that

5

where N₁ and N₂ are the normal matrices of the adjusted horizontal (N,E) coordinates at epochs t₁ and t₂, propagated from the GPS geocentric coordinates covariance matrices; n₁ and n₂ are the number of observations at epoch t₁ and t₂, respectively; m is the number of estimated parameters; and are the solution variances.

Eq. (5) holds under the hypothesis

6

which can be verified using the following test

7

In our case, this condition is verified for the 2000, 2001 and 2003 estimated variances at a significance level α = 10 per cent.

Eq. (5) is derived from the standard Fisher's test on least squares adjusted parameters that is commonly used in control problems (Mikhail 1976). It can be applied to the whole displacement vector (the entire network), on a part of it (e.g. one of the two transects) or on a single point. Moreover, it has been generalized to allow distinct normal matrices at epochs t₁ and t₂ and to be applied to horizontal components only (see also Anzidei et al. 1996). In fact, in the standard Fisher's test, the observation scheme and the stochastic model, that is, the design matrices and the observation covariance matrices, are supposed to be the same at the two control epochs. This assumption is not realistic in case of GPS measurements, since at least the satellite configuration is varying in time. Besides, as we are mainly interested in the horizontal displacements coming from the GPS network, we set up the test for the horizontal displacements , transforming the adjusted geocentric coordinates (X, Y, Z) to local coordinates (N, E, U) and propagating their covariances.

Under the hypothesis (6), a second statistic can also be used in order to rank the different model predictions versus the GPS estimated displacements.

The quadratic form of the normalized horizontal residuals U is chi-square distributed

8

where U is the difference between the horizontal GPS displacements, and the horizontal geophysical model displacements is the inverse covariance matrix of the horizontal displacement (N,E), propagated from the geocentric coordinate covariance matrix. The lower is the χ²_m value the better is the agreement with the tested model.

In order to verify the statistical significance of the geodetic measurements, the Fisher's test is here applied to the null-displacement model . The hypothesis

9

has been tested according to eq. (5). At a significance level α = 10 per cent, this hypothesis has been rejected both for the global network and the two transects separately, as shown in Table 3 where the values F_sp are much larger than the values F_th. Thus, the GPS derived horizontal displacements are statistically different from zero and it is meaningful to compare them with the predictions coming from the geophysical models.

Table 3.

Statistical significance of the GPS horizontal displacements. The theoretical Fisher value F_th has been computed at significance level α = 10 per cent. For each couple of campaigns, the first column stands for the whole network (All), the second column for the transect through the main fault (Tr.A, SPEL-SEFR) and the third column for the second transect (Tr.B, RASI-CENT).

To be used as reference for the χ²_m values listed in Sections 5 and 7, in Table 4 we list the χ² values for the null-displacement model, where we have dropped the subscript m in the notation.

Table 4.

χ² values for the null-displacement model.

5 Viscoelastic Relaxation

The main purpose of this section is to test the effect of different fault and earth models on the observed post-seismic deformation through viscoelastic relaxation only. For this purpose, the three different fault models (Table 1) are coupled to a number of candidate earth models (Table 5).

Table 5.

Crustal earth models used in this study. UC: upper crust, TZ: crustal transition zone, LC: lower crust.

Viscoelastic relaxation modelling is computed on the basis of upgraded normal mode relaxation models with a Maxwell viscoelastic rheology for a vertically stratified spherical Earth (Sabadini & Vermeersen 1997; Riva & Vermeersen 2002).

Though aware of the importance of power-law rheologies in controlling viscoelastic relaxation (e.g. Freed & Burgmann 2004), we decided to model relaxation by means of a linear rheology due to the characteristics of our data set:

• (1)

  Sparse campaign data: the density of the network (only 10 sites) and the availability of yearly campaigns does not allow to follow in detail the time-dependence of the relaxation process.

• (2)

  Specific time-window covered by the observations: we have at our disposal measurements between 2 and 5.5 yr after the earthquakes, which do not cover either the early post-seismic part (first 2 yr) or the relaxation occurring over time scales longer than 5 yr.

Moreover, the moderate size of the earthquakes can only lead to relaxation rates in the order of a few millimetres per year, so that measurement errors would partially mask qualitative differences in the relaxation process. Therefore, our findings cannot exclude non-linear or transient rheologies, although our viscosity values can be interpreted as indicative of those attained by power-law rheologies at various depths, in case non-linear flow is operative in the crust.

In Fig. 3, we show GPS displacements compared to one example of viscoelastic relaxation for years 2003–2000, for the three different fault models Zollo FM, Zollo-Basili FM and Salvi FM, with the candidate earth model TZ18.

Figure 3.

2003–2000 GPS and viscoelastic relaxation for model TZ18. Error ellipses represent the 98 per cent confidence level.

As far as the first transect, Tr.A (Fig. 1), is concerned, Zollo FM provides a definitely better fit than the geodetic Salvi FM, especially at the furthermost sites (SPEL, VALL and SEFR). The addition of the shallow nearly vertical fault extension by Basili & Meghraoui (2001) leads to a minor change in the fit, mostly effective at site MONT. All fault models predict a similar motion at CERE, whereas Salvi FM represents the best fit at COLL, even if all predictions fall within the error ellipse.

For the second transect, Tr.B (Fig. 1), the situation is different, because Zollo FM and Zollo-Basili FM better reproduce the two northernmost sites DIGN and CENT, whereas Salvi FM provides a better fit at POPO and RASI.

The general features of earth model TZ18 are maintained by the other models listed in Table 5. For this reason, we will discuss the performance of the various Earth and fault models on the basis of the Chi-square values listed in Table 6, as the hypothesis is not verified for any viscoelastic relaxation model, both globally and per transect.

Table 6.

χ² test for various Earth and fault models, for GPS campaigns 2003–2000.

Moreover, since the trend of all viscoelastic relaxation scenarios remains rather constant in the time-span observed by our GPS measurements, we have decided to compare viscoelastic models only against the campaigns 2003–2000. We have thus neglected the particular year-to-year behaviour mainly marked by a clear twisting in the vectors of the displacement field that cannot be reproduced by any viscoelastic relaxation model. However, we have verified that our χ² analysis, listed in Table 6, is also representative of the campaigns 2001–2000 and 2003–2001 when considered separately.

5.1 Sensitivity to TZ and LC viscosity

All listed earth models provide, for the seismological fault model Zollo FM, an improvement with respect to the null-displacement case (Table 4). In particular, a change in the LC viscosity has only a small impact on the expected viscoelastic relaxation: almost no difference is seen between a LC that is either elastic (mod.TZ18) or with a viscosity of 10¹⁸ Pa s (mod.LC18), whereas a further lowering of LC viscosity to 10¹⁷ Pa s (mod.LC17) leads to a worsening of the fit. According to this result, in most earth models we have considered an elastic LC, in order to isolate the effect of viscoelastic relaxation in the TZ, which appears to be the main contributor to the observed surface deformation. A change in TZ viscosity from the value of 10¹⁸ Pa s (mod. TZ517 and TZ518) leads to a gradual worsening of the fit, due to a general reduction in the magnitude of viscoelastic relaxation. In fact, a higher viscosity slows down the relaxation process, whereas a lower viscosity shifts most of the displacement to the first 3 yr after the earthquakes, thus prior to the considered GPS campaigns.

The addition of a shallow and high-angle fault, as proposed in Basili & Meghraoui (2001), produces a similar behaviour with respect to changes in the earth models, but the weight of the two fault segments changes. In fact, a general deterioration of the fit for the first transect is accompanied by an improvement in the fit of the second transect; nonetheless, the global fit is slightly worse than that obtained with Zollo FM.

The situation is different for the geodetic model Salvi FM that points to a TZ viscosity of 5 × 10¹⁷ Pa s. However, the agreement of this fault model with the GPS measurements is generally worse than that obtained with the seismological model, with χ² values mostly close to the results of the null-displacement model (Table 4).

5.2 Sensitivity to TZ thickness and location

An important issue about the earth models is represented by the actual thickness and location of the TZ. From Section 5.1, we learned that the most important contribution to the viscoelastic relaxation is coming from the low-viscosity TZ, which has been so far represented by a 12-km-thick layer underlying an 8-km-thick elastic upper crust (see Table 5). An important issue about the earth models is represented by the actual thickness and location of the TZ.

The hypothesis of an 8-km-thick UC is justified by a number of arguments:

• (1)

  The abrupt cut-off of the 1997 aftershock sequence (e.g. Amato et al. 1998; Cattaneo et al. 2000).

• (2)

  The shallow regional seismicity (Chiarabba et al. 2005).

• (3)

  Results from deep crust reflection studies (Pialli et al. 1998) and more recent shear wave velocity models (Chimera et al. 2003).

Nonetheless, it is worth analysing the effect of a slightly thicker UC, as in model ‘d10–d20’: the result for the seismological fault model is a marked reduction in the magnitude of deformation, with the consequent deterioration of the global fit, whereas in the case of the geodetic fault model it leads to an improvement of the fit for the first transect and a deterioration for the second transect, leaving an almost unchanged overall fit.

On the other side, the location of the boundary between the TZ and the LC, so far fixed at 20 km depth, can reasonably be as shallow as 13 km, according to what is suggested by local seismic profiles Pialli et al. (1998) and S-wave velocity profiles Chimera et al. (2003). The earth model with a 5-km-thick TZ is labelled ‘d8–d13’: it provides a small increase of the fit for Salvi FM, mainly due to a better reproduction of the motion along the second transect, and a little deterioration in the fit of the first transect for Zollo FM, while the result of Zollo-Basili FM remains almost unvaried.

We can thus state that the thickness of the top elastic layer is a crucial parameter to reproduce the observed motions by means of viscoelastic relaxation, whereas the model is rather insensitive to the actual thickness of the first viscoelastic layer.

In conclusion, the earth model for viscoelastic relaxation that provides the best fit to the GPS measurements presents a low viscosity layer (TZ, η = 10¹⁸ Pa s) located below a rather thin UC (8 km thick). The thickness of the TZ (between 5 and 12 km) and the viscosity of the LC have only a small impact on the modelled deformation and are thus not sufficiently resolved by the data. Moreover, the preferred fault model is represented by the solution published by Zollo et al. (1999).

6 Poro-Elastic Relaxation

Besides viscoelastic relaxation, another potentially important Post-seismic process is represented by poro-elastic relaxation, due to changes in pore-fluid pressure induced by the earthquakes (e.g. Peltzer et al. 1998; Fialko 2004). We computed the fully relaxed poro-elastic signal by calculating the difference between the elastic deformation induced by the earthquakes in undrained and drained crustal rocks. In the attempt of maximizing the poro-elastic response, we have chosen extreme values of Poisson's ratio, namely ν = 0.35 and ν = 0.20, to represent the undrained and drained elastic moduli, respectively. The elastic deformation has been modelled by means of an Okada (1985) half-space model (Feigl & Dupre 1998).

In Fig. 4, we show the resulting displacements vectors at the GPS sites, together with the observed displacement for years 2003–2000, for the case of Zollo FM. As we can see, the signal induced by poro-elastic relaxation is much smaller than the observed GPS motion, and this result becomes even more clear when we consider two important aspects:

Figure 4.

2003–2000 GPS (black arrows) and total poro-elastic relaxation (red arrows). Error ellipses represent the 98 per cent confidence level.

• (1)

  In the first place, the chosen Poisson's ratios represent extreme values, so that actual poro-elastic relaxation can easily be smaller than what we have computed.

• (2)

  Second and most important, our GPS motions concern a specific time-window, between 3 and 5.5 yr after the earthquakes, so that we are not able to observe the large portion of poro-elastic relaxation that has probably taken place right after the earthquakes, as also inferred from the rock properties of the UC derived from 3-D tomography (Monna et al. 2003).

Considering those results, in the rest of our study we have decided to regard poro-elastic relaxation as a second order process; therefore, we have not included it further in our effort of modelling the observed GPS motions.

7 After-Slip

We have discussed how viscoelastic and poro-elastic relaxation are largely underestimating the horizontal GPS displacements. Moreover, the measured displacements of the near fault sites, shown in Fig. 2, present large variations in the motion directions between campaigns 2001–2000 and 2003–2001 that are not reproduced by any relaxation model. In order to account for the missing deformation, in this section we explore the contribution of after-slip, another potentially important post-seismic process (e.g. Savage et al. 1994; Pollitz et al. 1998). By after-slip we mean potentially either shallow fault creep, aseismic slip below the fault zone, or both.

We have studied the effect of after-slip on separate segments for different combinations of the two fault models (Zollo FM, Salvi FM) for the two main shocks. The after-slip segments are as follows.

• (1)

  Coseismic faults.

• (2)

  Updip extensions of the faults, with the same dip.

• (3)

  Vertical updip extension of the faults, as proposed by Basili & Meghraoui (2001).

• (4)

  Downdip extension of the faults.

• (5)

  Horizontal faults at the base of the UC.

For each fault segment, we explore variable amount of slip, rake and fault width. By means of a trial-and-error approach, minimizing the misfit between measurements and model results, we search for the best values for the above parameters.

Slip on the coseismic faults or on their downdip extentions does not provide any satisfactory fit to the observed GPS motions. Our tests have shown how the effect of accelerated creep on and below the coseismic faults is only relevant for near-fault sites on the hanging wall, but also in this case the impact is reduced to negligible values by our trial and error minimization of the misfit with the measurements.

Differently, slip on shallow fault extensions, besides being consistent with studies on the mechanics of after-slip (Marone et al. 1991), provides an important contribution to the observed motions, in particular to the reversal of the near-fault displacement observed after 2001, and will be discussed extensively.

We start our analysis with the displacements in the years 2001–2000 for Zollo FM, shown in Fig. 5. It is evident from panel (a), where GPS motions are compared with the predictions of viscoelastic relaxation for earth model TZ18, that most of the deformation remains unaccounted for, especially at the three sites closer to the main fault (CERE, COLL and MONT). Large near-fault motions are only reproduced by a rather shallow slip. Moreover, the direction of motion requires that slip occurs on a low angle fault, likely corresponding with an updip extension of the main shock rupture. The best-fitting model of shallow after-slip, shown in panel (c), allocates about 7 cm normal slip with a left-lateral component on the updip extensions of both faults, between a depth of 2 and 4 km. A better approximation to the motions observed at the furthermost sites requires the addition normal slip on a horizontal fault plane at the base of the UC. Keeping fixed the length of the horizontal segment, equal to the length of the coseismic fault, we have tested different locations, from the hanging wall to the footwall, and along-dip dimensions, from 3 to 18 km: the best-fitting model, shown in panel (b), presents 5 cm of normal slip on a 6-km-wide fault located below the footwall of the main fault.

Figure 5.

2001–2000 GPS (black arrows) and predicted displacements (coloured arrows) for Zollo FM. Error ellipses represent the 98 per cent confidence level. Panel (a) represents viscoelastic relaxation for earth model TZ18; panel (b) after-slip on a 6-km-wide horizontal fault at the base of the UC, under the footwall of the main fault; panel (c) after-slip on the updip extension of both faults, between 2 and 4 km depth; panel (d) the summation of the three contributions from panels (a), (b) and (c).

The three contributions, namely viscoelastic relaxation, slip on the updip extension and deep horizontal after-slip, are summed to provide the best-fitting model, represented in panel (d) of Fig. 5. We clearly see that a good agreement between measured and modelled displacements is reached at most sites, with the exception of those sites that show a deviation from the general trend of motion. The fact that we apply a homogeneous slip on the various fault segments is likely the main reason for the trade-off between the different contributions, that prevents us from obtaining a better fit at some specific sites without a general deterioration of the fit for the rest of the network. The improvement with respect to pure viscoelastic relaxation is also demonstrated by the χ² values listed in Table 7.

Table 7.

χ² test for the combined viscoelastic relaxation—after-slip models. Sites passing the Fisher test at significance level α = 1 per cent are also listed.

The effect of after-slip for Salvi FM is shown in Fig. 6. Differently from the previous case, the most important contribution comes from 9 cm of normal slip on a 9-km-wide horizontal fault plane at the base of the UC, displayed in panel (b). This segment has been localized below the footwall, according to the same procedure previously discussed for Zollo FM. A further refinement comes from the addition of 4 cm of normal slip on the upper segment of the main fault, between a depth of 1.5 and 0.5 km, displayed in panel (c).

Figure 6.

2001–2000 GPS (black arrows) and predicted displacements (coloured arrows) for Salvi FM. Panels (a) and (d): same as Fig. 5. Panel (b): after-slip on a 9-km-wide horizontal fault at the base of the UC, under the footwall of the main fault; panel (c): after-slip on the upper segment of the main fault, between 0.5 and 1.5 km depth.

The χ² values listed in Table 7 show how the best model for Salvi FM has a slightly worse overall fit that Zollo FM. In particular, Zollo FM has a better performance at COLL, POPO and DIGN, whereas Salvi FM provides better results at RASI and CENT.

GPS motion vectors for the years 2003–2001, as already anticipated, present rather large differences with respect to years 2001–2000: CERE, COLL and MONT invert the motion direction, while POPO, DIGN and CENT are rotated by about 90° clockwise. Vectors are displayed in Fig. 7, together with model results for Zollo FM.

Figure 7.

2003–2001 GPS (black arrows) and predicted displacements (coloured arrows) for Zollo FM. Error ellipses represent the 98 per cent confidence level. Panel (a) represents viscoelastic relaxation for earth model TZ18; panel (b) after-slip on a 4-km-wide horizontal fault at the base of the UC, under the hanging wall of both faults; panel (c) after-slip on a nearly vertical updip extension the main fault, between 2 and 4 km depth; panel (d) the summation of the three contributions from panels (a), (b) and (c).

The motion of sites CERE, COLL and MONT can be grossly explained by 10 cm of after-slip on the almost vertical updip extension proposed by Basili & Meghraoui (2001), here confined between a depth of 4 and 2 km on the main fault, and shown in panel (c). Slip on a horizontal fault at the basis of the UC becomes more important, with the allocation of about 10 cm of normal slip on a 4 km-wide fault below the hanging wall for both faults. It provides necessary motion at sites SPEL, VALL, DIGN and CENT, and at the same time contrasts the otherwise exceedingly large NE-motions at sites MONT and SEFR, as shown in panel (b). Magnitudes and directions at SPEL, VALL and SEFR are similar to those of viscoelastic relaxation, represented in panel (a). The summation of the three contributions, displayed in panel (d) of Fig. 7, represents our best-fitting model: part of the deformation remains unexplained, particularly at MONT, where the two vectors are almost perpendicular, and POPO, where the large measured displacement is completely unaccounted for. Those two sites, however, are at the border of an abrupt change in motion, namely the observed shortening of the baselines MONT-SEFR and RASI-POPO, which represent a small scale behaviour difficult to reproduce with our approach, where homogeneous slip is applied on relatively large fault segments. Again, the significant improvement with respect to pure viscoelastic relaxation is seen in the χ² values listed in Table 7.

It proves difficult to construct an adequate after-slip model for 2003–2001 starting from the coseismic rupture geometry of Salvi FM. In this case, in fact, we miss the possibility of allocating slip on a shallow and nearly vertical fault, since the results published by Basili & Meghraoui (2001) have been obtained after using as reference rupture model the deeper fault of Zollo et al. (1999). From the χ² values listed in Table 7, we can see how only a minor improvement to the fit of viscoelastic relaxation comes from the addition of the best after-slip model for Salvi FM, represented by 3 cm of normal slip on the main fault between a depth of 3.5 and 1.5 km.

The rather good agreement between the GPS motions and the predictions of the best post-seismic deformation models allows to go beyond the χ² values, which are only a measure of the misfit, and discuss the statistical significance of the model predictions by means of the Fisher's test described in Section 4.3. In the right column of Table 7, we list the sites for which the model predictions pass the Fisher's test at significance level α = 1 per cent.

In the case of the results for Zollo FM in years 2001–2000, we see that viscoelastic relaxation alone provides a significant fit only at sites SEFR, SPEL and DIGN. After-slip alone allows successful predictions at five sites on Tr.A, missing CERE, and the same result is obtained by the combined after-slip and relaxation model, meaning that the two models are statistically equivalent. Larger differences are present on Tr.B, because both models fit DIGN, but after-slip alone fits also POPO and CENT, whereas the combined after-slip and relaxation model fits RASI.

Salvi FM for years 2001–2000 provides similar results, further missing only COLL on Tr.A and DIGN on Tr.B. Both sites are situated directly above the two main faults according to Salvi et al. (2000): therefore, the observed misfit might be due to an erroneous fault location either horizontally or in depth.

For the 2003–2001 campaigns only Zollo FM provides a statistically significant fit to the GPS data. Pure relaxation fits three sites on Tr.A, SEFR, MONT and SPEL and RASI on Tr.B. After-slip alone and the combined after-slip relaxation model are statistically equivalent and fit the four southernmost sites of Tr.A, in addition to RASI and DIGN on Tr.B. The degrading of the fit at MONT and SEFR introduced by the after-slip model is due to the necessity of fitting large and opposite motions at the remaining four sites of Tr.A and could not be avoided. On Tr.B, only DIGN shows a motion that is both significantly different from zero and reproduced by the after-slip model.

In Fig. 8, we show rates of horizontal baseline variations along Tr.A, with respect to the westernmost site SPEL. GPS results for campaigns 2000–1999, represented by open stars joined by a dotted line, have only a reference purpose because they have not been used in the rest of the study. Model results are obtained by the joint contribution of relaxation and after-slip for fault model Zollo FM. The fit between measurements and model results is in most cases well within the 68 per cent confidence level. In agreement with what has already been largely discussed in the previous sections, also in the baseline representation we can see how SPEL-VALL is always shortening, at a rate decreasing after 2001, while the three sites nearest to the fault are moving towards SPEL between 2000 and 2001 (green dots and squares) and in the opposite direction between 2001 and 2003 (red dots and diamonds). The baseline SPEL-SEFR is not showing any significant motion between 1999 and 2003, while model results give extention rates of 2–3 mm yr⁻¹, but this discrepancy is consistent with the expected measurement error. The large variation in the direction of motion between 2000 and 2003 can only be reproduced by after-slip, since viscoelastic relaxation alone gives rise to a deformation pattern analogous to the one between 1999 and 2000, as reported by Aoudia et al. (2003).

Figure 8.

Horizontal baseline variation rates along Tr.A. Model displacements represent the best-fitting result of combined relaxation and after-slip for Zollo FM. Error bars represent the 68 per cent confidence level.

8 Discussion

The study of deformation between 1999 and 2003 has clearly indicated that long-lasting after-slip is the main post-seismic process responsible of the observed GPS motions, in terms of both magnitudes and directions. This conclusion is made possible by the extension of the GPS network from year 2000, which allows to obtain reliable planar displacements, as testified by the positive Fisher's test values listed in Table 3.

The best results are obtained allowing various patches adjacent to the two main faults to slip aseismically, using Zollo FM as preferred coseismic model. In particular, we need shallow slip to match the near fault sites and slip at depth to fit the far-fault sites.

The observed post-seismic deformation requires the contribution of both faults, although the fault reactivated by the second large event (9:40, M_w 6.0) is the leading one. Regarding this later fault, the main preferred coseismic model would require a rupture depth as in Zollo FM that is deeper than Salvi FM. Moreover, the fact that the preferred fault model ruptures the lower half of the seismogenic UC is a key element to justify the activation of the nearly vertical updip extension initially proposed by Basili & Meghraoui (2001), and the occurrence of slip at the base of the UC itself. For the fault reactivated by the first event (0:33, Mw 5.7), the main difference between Zollo FM and Salvi FM stands in the location of the fault projection at the surface, where Zollo FM is at about 4 km to the NE with respect to Salvi FM. The fit between the observations and model predictions for the pertinent transect as a whole (Tr.B) does not allow us to choose between Salvi FM and Zollo FM. However, the motion of the near fault sites (e.g. DIGN, Fig. 5d versus Fig. 6d) gives more weight to Zollo FM.

Between 2000 and 2001, the shallow component is located on the up-dip extension of both faults between a depth of 2 and 4 km and accommodates about 7 cm of normal slip. At depth, we allow the base of the footwall of the main fault to slip by 5 cm in normal direction above the low-viscosity TZ. The required equivalent moment for this 1-yr period amounts to at least 12 per cent of the total coseismic moment of the two events.

Between 2001 and 2003, the very different deformation pattern for the near fault sites requires a considerable amount of normal slip, about 10 cm, to be located on a nearly vertical updip extension of the main fault, between a depth of 2 and 4 km. At depth, the activated fault plane shifts SW under the hanging-wall for both faults, allocating 10 cm of normal slip. The required equivalent moment release for this 1.6-yr period amounts to at least 15 per cent of the coseismic moment. Therefore, the rates of after-slip on the faults are slightly decreasing from 2000 to 2003.

The last component of after-slip, allocated on a horizontal fault at the base of the seismogenic layer, could be considered physically analogous to viscous relaxation below the brittle-ductile transition. However, a comparison between this contribution and the best model of viscoelastic relaxation [panels (a) and (b) of Figs 5–7] shows how the two processes affect the motion of the GPS sites in different ways. Therefore, in the specific case, we do not regard the two processes as being equivalent.

In spite of the important role of after-slip, the long-wavelength deformation exhibited by the far-fault sites requires a contribution of viscoelastic relaxation. The best results for viscoelastic relaxation are obtained with Zollo FM, when a layer with viscosity of 10¹⁸ Pa s (TZ) is located below an 8-km-thick upper crust. The observation window of our GPS measurements spans between 2 and 5.6 yr after the earthquakes. A much longer observation time would help to detect the presence of any higher viscosity, while the effect of a much lower viscosity has probably extinguished in the first year after the earthquakes.

We have verified that another potentially active post-seismic deformation process, namely poro-elastic relaxation, is not capable of affecting our GPS sites motions significantly.

9 Conclusions

We have shown that the campaign occupation of a small size GPS network is capable of detecting deformation signals of the order of a few millimetres per year, provided that accurate processing and antenna positioning is realized.

The comparison between GPS measurements and displacement predictions coming from different models of post-seismic deformation allows to put some constraints on the earth structure, the rupture models and the specific post-seismic processes active in the area of the 1997–1998 Umbria-Marche earthquake sequence.

The fit to the GPS motions is obtained when using the coseismic fault model published by Zollo et al. (1999) and allowing various patches adjacent to the two main faults to slip aseismically. The slip on different shallow patches reflects an important on-going fault creep, while the after-slip, at the base of the seismogenic UC, reflects the presence of a possible discontinuity at 8 km. This discontinuity is likely rheological, decoupling the seismogenic UC from the TZ, and is the locus of a basal shear favouring localized time-dependent deformation. The preferred rheological model, obtained from viscoelastic relaxation modelling, consist of an 8-km-thick elastic UC, underlined by a 12-km-thick TZ with a viscosity of 10¹⁸ Pa s. Contributions from deeper layers, LC and upper mantle, are negligible for the time-span of the observations.

The first order process governing the observed post-seismic deformation following the Umbria-Marche earthquakes is controlled by after-slip both above and below the fault zones.

Acknowledgments

This work is fully supported by the Italian MIUR-PRIN 2004 project: Active deformation at the northern boundary of Adria. We thank the staff of Politecnico of Milan, University of Milan and University of Trieste for the help during the GPS campaigns. We would also like to thank M. Crespi and M. Meghraoui for their comments on this paper.

Author notes

*

Now at: Geophysics Section, Department of Earth Sciences ‘A. Desio’, University of Milan, Italy.