Rupture history of the 1997 umbria-Marche (central Italy) main shocks from the inversion of GPS, DInSAR and near field strong motion data

We investigate the rupture history of the three largest magnitude earthquakes of the 1997 Umbria-Marche sequence by inverting GPS, DInSAR and near-source strong motion waveforms. We use the frequency domain in- version procedure proposed by Cotton and Campillo (1995) and calculate the Green’s functions for a layered half-space using the discrete wavenumber and reflectivity methods. We first invert GPS measurements and DInSAR interferograms to image the coseismic slip distribution on the fault planes in a layered half space for the two earthquakes that occurred on September 26, 1997 at 00:33 UTC ( M w = 5.7) and 09:40 UTC ( M w = 6.0) near Colfiori- to. We also invert DInSAR interferograms to infer the slip distribution during the subsequent earthquake that occurred on October 14, 1997 at 15:23 UTC ( M w = 5.6) in the SE section of the seismogenic zone near Sellano. We also explore the set of acceptable solutions using a genetic algorithm to have information on the available resolution of geodetic data. The slip models obtained by geodetic data inversion are used to perform a forward mod- eling of strong motion waveforms for all three events. We adopt a constant rupture velocity of 2.6 km/s and a constant rise time of 1 s. Our results show that these rupture models provide an acceptable fit to recorded waveforms. Finally, we invert the recorded ground displacements, collected during the September 26th 09:40 main shock and the October 14th Sellano earthquake, to constrain the rupture history. We use the geodetic slip distribution as start- ing model for the iterative inversion procedure. The retrieved rupture models are consistent with those inferred from geodetic data and yield a good fit to recorded seismograms. These rupture models are characterized by a het- erogeneous slip distribution and an evident rupture directivity in agreement with previous observations. Abstract The scaling law of the seismic spectrum experimentally calculated at Mt. Vesuvius and Campi Flegrei is used to con- strain the estimate of the maximum expected peak acceleration of ground motion. The scaling law was calculated for earthquakes recorded at BKE and OVO stations in the period 1997-2000 for Mt. Vesuvius and for earthquakes oc- curred during the 1983-1984 bradyseismic crises at Campi Flegrei. For Mt. Vesuvius the scaling law clearly deviates from a constant stress drop relation in the whole range of magnitude (0.4 < M D < 3.6) whilst constant stress-drop is found for Campi Flegrei data (0.7 < M D < 3.4). These results are used to give a first estimate of the maximum ground motion corresponding to the largest magnitude ( M max ) inferred in the two investigated areas by the Gutenberg-Richter formula. The values of the seismic moment M 0 and the characteristic source radius corresponding to M max are used to evaluate the peak ground acceleration PGA. This parameter is determined by stochastic simulation of ground motion. Two different methods (Random Vibration Theory (RVT) (Boore, 2003) and ground motion generated from a Gauss-ian distribution with σ = a rms (GMG)) give slightly different values of PGA. The values of PGA were 0.10 g (RVT), 0.14 g (GMG) for Mt. Vesuvius and 0.04 g (RVT and GMG) for Campi Flegrei.


Introduction
The 1997-1998 Umbria-Marche (central Italy) seismic sequence struck the northern part of the Apennines and it is characterized by six earthquakes with magnitude larger than 5.0 (see Amato et al., 1998;Ekström et al., 1998) and more than 2000 smaller shocks (Deschamps et al., 2000;Chiaraluce et al., 2003). Most of the total seismic moment was released within the first two months of activity , but seismicity lasted until the end of 1998 (Deschamps et al., 2000). The seismic sequence began with a M L 4.3 foreshock on September 3rd 1997, which was followed by several aftershocks within the subsequent two weeks (see Ripepe et al., 2000). The largest shocks of the sequence occurred on September 26 (two events at 00:33 and 09:40 UTC) and October 14 (15:23 UTC) and caused severe damage and several causalities. These seismic events allowed the recording of the largest peak ground accelerations (nearly 0.6 g at the Nocera Umbra and Colfiorito sites, see Bouchon et al., 2000;Cultrera et al., 2003) within the Italian strong motion data set. The evolution of seismicity and earthquake locations are described by Deschamps et al. (2000) and Chiaraluce et al. (2003), while the details of the fault geometries and the faulting mechanisms are discussed in Chiaraluce et al. (2004). Ekstrom et al. (1998) calculated CMT fault plane solutions of the largest shocks of the sequence, while Chiaraluce et al. (2003) calculated the focal mechanisms of aftershocks using relocated hypocenters through a 3D local tomographic inversion of travel times (Chiarabba and Amato, 2003) and first motion polarity data.
On September 26, 1997, two closely located earthquakes (see fig. 1) struck the Colfiorito area within nine hours from each other rupturing in opposite directions two normal faults (Amato et al., 1998). The first one (Mw 5.7) occurred at 00:33 UTC and broke toward the SE a SW dipping normal fault, while the second one (Mw 6.0) occurred at 09:40 UTC and ruptured toward the NW a distinct fault segment dipping to the SW (see Chiaraluce et al., 2003Chiaraluce et al., , 2004. This second shock caused severe damage to the historical heritage of the area and in particular it caused the collapse of part of Assisi Cathedral. The seismic activity continued after these two main shocks and seismicity activated a series of normal faults oriented along a NW-SE direction along the Apennines (see fig. 1) and dipping to the SW (see Chiaraluce et al., 2003Chiaraluce et al., , 2004. On October 14, 1997, another major seismic event (M w 5.6) occurred in the SE section of the seismogenic zone, near Sellano (event 7 in fig. 1). The fault plane solutions of the earthquakes that occurred during the 1997-1998 seismic sequence are consistent with the active horizontal extension perpendicular to the Apenninic chain (Frepoli and Amato, 1987;Montone et al., 1997;Chiaraluce et al., 2004).
The slip distribution during the three largest shocks of the sequence was previously investigated through a forward modelling of GPS measurements by Hunstad et al. (1999) and by GPS and DInSAR data by Salvi et al. (2000). The fit to geodetic data allowed the fault dimensions and geometry to be constrained. More recently, Lundgren and Stramondo (2002) imaged the slip distribution for the September 26 earthquake by inverting GPS and DInSAR data using the simulated annealing numerical procedure. These authors used a single fault plane for the two shocks of September 26 and concluded that slip is concentrated at depth ranging between 4 and 6 km, with a smooth gradient toward the surface. Regional broad-band seismograms were analysed by Pino and Mazza (1999) to constrain rupture directivity for the largest shocks: this study confirmed that the two shocks of September 26 ruptured two distinct normal faults in opposite directions. Strong motion data were modelled by Zollo et al. (1999) and Capuano et al. (2000) to image the gross feature of the rupture model for the September 26 (09:40) main shock. All these investigations provided consistent results. However, the coseismic GPS displacements, the surface deformation pattern resulting from DInSAR observations and the radiated seismograms have only been modelled and interpreted separately.
The aim of this paper is to study the rupture history of the three main shocks of the 1997-1998 sequence by the inversion of geodetic (GPS and DInSAR) data and strong motion seismograms. Despite their moderate magnitude, these events allowed the collection of a multidisciplinary and accurate data set, which are very useful to constrain the slip distribution and the rupture propagation on the three main fault segments. The study of such moderate size events is relevant for different reasons. From a methodological point of view, it is important to know the possibility of studying the rupture process of moderate-magnitude events, rupturing fault planes of less than 100 km 2 , using near field data. For the assessment of seismic haz-ard, this source study is important to understand the damage distribution and to discriminate the effects ascribed to local site geology or to the source itself. It is also important, in the case of a sequence composed of several earthquakes with similar magnitude, to know the details of fault geometry and slip distribution to attempt to understand the interactions between distinct events. Cocco et al. (2000), for instance, investigated fault interaction through elastic stress transfer among the mainshocks of the Colfiorito sequence adopting a homogeneous distribution of slip on the fault planes. The study of the rupture history is also of relevance for interpreting the regional tectonics, the fault locations and the surface observations (see Cinti et al., 1999Cinti et al., , 2000Chiaraluce et al., 2004 and references therein).
Most of the damaging earthquakes that occurred recently in Europe (1980 Irpinia, Italy;1983 Corinth and Athens, Greece) are normal faulting events. Studies aimed at the inversion of seismological and geodetic data for such events are less frequent than those concerning strike slip events (Imperial Valley, 1979;Landers, 1989;Loma Prieta, 1992;Hector Mine, 1999;Izmit, 1999;Tottori, 2000). The study of the rupture process and slip distribution for normal faulting earthquakes is fundamental for seismic hazard in Italy.

The available data
In this study we invert geodetic and strong motion data to constrain a kinematic rupture model for the largest events of the Umbria-Marche seismic sequence. The geodetic data set comprises the deformation pattern resulting from Differential Interferometry Synthetic Aperture Radar (DInSAR) and the vertical and horizontal displacements measured through the Global Positioning System (GPS). To image the rupture propagation on the fault plane we use near field strong ground motions recorded close to the faults and collected by different Italian institutions. Fig. 1. Map of seismicity during the 1997 Colfiorito sequence. The epicenters and focal mechanisms of the main shocks (shown in lower hemisphere projection) are indicated with solid circles. The source process of the 3 strongest earthquakes of the sequence (namely, the Colfiorito event # 2, occurred on September 26, 1997, at 00:33 UTC; the Colfiorito event # 3, also occurred on September 26 at 09:40 UTC and the Sellano event # 7, occurred on October 14 15:23 UTC) are investigated in this study by using GPS, DInSAR and strong ground motion observations.

DInSAR interferograms
In this paper we use ERS1-ERS2 SAR images. Raw data were processed by INGV and CNR-IRECE research groups to produce coseismic interferograms relative to one or several earthquakes of the sequence (Stramondo et al., 1999). The resulting interferograms provide a fringe pattern (see fig. 2a,b) that describes the surface displacement field projected onto the slant-range direction (i.e. the satellite-to-ground line of sight). Each fringe corresponds to a displacement of 28.3 mm in slant range. The various applications of radar interferometry to changes in Earth's surface are described by Massonnet and Feigl (1998). After the first spectacular results relative to the Landers earthquake (Massonnet et al., 1993), the DInSAR technique has been demonstrated to be a powerful tool to study surface deformations due to earthquakes. The first interferogram obtained during the Umbria-Marche se-quence is described and interpreted by Stramondo et al. (1999); Salvi et al. (2000) analyzed new ascending and descending interferograms to improve the interpretation of SAR data. Lundgren and Stramondo (2002) completed the analysis of DInSAR data with a new interferogram showing the surface displacement pattern caused by the October 14th earthquake.
In the present paper we limit our study to the interpretation of 2 interferograms. The first one is obtained by a 35-day ERS2 descending image pair (Stramondo et al., 1999) and covers the period between September 7 and October 12. The altitude of ambiguity (ha, i.e., the altitude variation corresponding to a phase change of 2π) is about 65 m. Fringes observed in this interferogram are mainly due to the deformation field of the 2 main events of September 26 occurred at 00:33 and 09:40 UTC. We digitized fringes obtaining 194 range change values each one every 0.5 km ( fig. 2a). The extracted points reveal an Fig. 2a,b. a) Fringe pattern extracted from available interferograms (see text) for the Colfiorito area. Each fringe is a contour line of equal displacement of the ground surface in the slant range direction, representing a displacement increment of 28.3 mm measured between September 7, 1997 and October 12, 1997. b) Fringe pattern extracted from the interferograms for the Sellano area. It represents the displacement field in the ground to satellite direction observed between August 9, 1997 and October 17, 1997. increase in the distance between the earth surface and the satellite. This observation corresponds to a deflation of the ground surface above the fault consistent with the normal faulting CMT solutions and the extensional tectonics of the region. Unfortunately, no SAR images are available for the period comprised between the two main events of September 26. This makes it necessary to consider the two events in a single model. We finally use a second interferogram representing the deformation field in the ground to satellite axis between August 9 and October 17 (Salvi et al., 2000). It is a ERS1-ERS2 70day ascending interferogram with ha 65 m. We used the 300 values extracted from fringes and inverted these data points relative to the Sellano region to access to the slip distribution on the fault due to the October 14 event ( fig. 2b).

GPS coseismic displacements
Several GPS benchmarks were set up and measured between 1992 and 1996 in Italy by the Istituto Geografico Militare (IGMI) and the stations located in the epicentral region were measured again in 1997 by the Istituto Nazionale di Geofisica (ING) during and after the sequence. The details of the GPS networks (IGM95 and TYRGEONET), the data processing and the computed coseismic displacement caused by the Umbria-Marche earthquakes are described and discussed by Anzidei et al. (2000) and modelled by Hunstad et al. (1999), Stramondo et al. (1999) and Salvi et al. (2000). Figure 3 shows the distribution of the GPS sites, used in previous studies to evaluate the dislocation model for the two earthquakes of September 26, the observed horizontal displacements and the associated errors (see Hunstad et al., 1999, and references therein). No relevant coseismic deformation was measured at the GPS sites for the October 14, Sellano, earthquake.

Strong ground motions
In this study we use the strong ground motion data recorded within 30 km from the epicenter by the accelerometers belonging to the Italian na-tional accelerometric network deployed by the ENEL (the Italian electric company), by ENEA and by the Servizio Sismico Nazionale (SSN). The stations that provided digital data are listed in table I, their distribution and the position of the fault planes that ruptured during distinct earthquakes are shown in fig. 4. Because of the difficulty in modeling high frequency radiation, the waveforms are integrated twice to get the ground displacement. The seismograms are then bandpass filtered between 0.1 and 1.5 Hz.
An important problem we had to face in modeling these data is the lack of absolute times of the recorded time histories. This information is crucial to perform an appropriate inversion for the rupture time in order to constrain the rupture front propagation on the fault plane. To estimate the absolute time for each time history we computed the theoretical arrival time of the direct P-and S-waves from the hypocenter to the recording sites, using the velocity structure listed in table II. We therefore picked the Pwave arrival times when they are present on the unfiltered accelerograms; they are listed in table III for each of the earthquakes investigated in this study. The error in absolute time is of the order of 0.2 s for digital stations, but it can be of the order of 1 s or larger for some analogical accelerograms. Thus, we have estimated the delay between triggering time from the earthquake origin time as a difference between the theoretical and the observed travel time.   The resulting triggering times for the selected ground motion time histories are listed in table III for each of the three largest shocks of the sequence. We analyze a time window including the P-and S-wave field for each waveform, whose length mainly depends on the distance between the fault and the considered station. The fit to the digital accelerograms recorded at Assisi (AS010) and Cerreto di Spoleto (CERT) is particularly important to investigate the rupture history at low frequencies because they are not affected by site amplification effects. We use other stations with lower weights, when the accelerograms were incomplete or perturbed by important site effects. This is particularly true for the Colfiorito recording site, where site effects dominate the ground motion time histories over a broad frequency band (see Rovelli et al., 2001;Di Giulio et al., 2003). The accelerograms recorded at Nocera Umbra are also affected by major site amplification effects, but their resonance alters ground motion amplitudes (see Cultrera et al., 2003 and references therein) at frequencies higher (f > 2 Hz) than those investigated in this study (f < 1.5 Hz).

Green's function calculation and crustal model
The Green functions are calculated for a layered half-space using the discrete wavenumber integration method (Bouchon, 1981) associated with the reflection transmission matrix method (Kennett, 1983). The crustal velocity model adopted in this study (see table II) consists of 3 layers over a half space, used by Deschamps et al. (2000) for aftershock locations. SAR and GPS observations. It is represented as a linear sum of n subfault contributions where Gik (f = 0) represents the static ground deformation (f = 0) for a unit constant slip on the subfault k (i.e., the static Green's Functions).

Inversion procedure
The inverse problem consists in minimizing the difference between the recorded data and those simulated through a forward modeling. The unknown parameters in our inverse approach are the slip amplitude, the rise time and the rupture time associated with each subfault. However, because we investigate moderate magnitude-events with relatively small fault planes, we do not invert for the rise time, which is taken constant during the inversion. We therefore invert geodetic and strong motion waveforms to image the slip and rupture time distributions on the fault plane.
The inversion of geodetic data assumes that the parameter vector m and the data vector d are linearly related by the operator G as d = Gm. The elements of d are the GPS measurements in the east, north and vertical directions and the DInSAR deformations in the ground to satellite axis. We follow Tarantola's (1987) formulation to calculate the model expectation m where Cm and Cd are the covariance matrix for m and d and m0 is the initial model. Seismological data are inverted by using the inversion approach proposed by Cotton and Campillo (1995), which is based on the non-linear inversion in frequency domain of strong motion recordings integrated to obtain ground displacements. We invert the complete Green's functions, therefore accounting for near field effects. The parameter vector m and the data We use this crustal model to compute both the static deformations (DInSAR and GPS observations) and synthetic seismograms in the near source for all the three largest earthquakes. The variation of the shear modulus with depth is also taken into account to compute the seismic moment.
In the literature, many inversions of geodetic data are usually performed in a homogeneous half space. However, the calculations performed considering a homogeneous half space yield significantly different results with respect to a stratified medium as shown by Savage (1998) and Cattin et al. (1999) for the deformation field and by Belardinelli et al. (1999) for the dynamic stress changes caused by the Irpinia earthquake. We use in this study the same layered crustal model to invert both geodetic data and ground motion waveforms.

Methodology
We calculate the synthetic ground motion waveforms in the frequency domain using the methodology proposed by Bouchon (1981) and Cotton and Campillo (1995). We express the Fourier transform of the seismogram Vi at a given (i-th) receiver as a linear sum of n subfault contributions where Gik(f) represents the Green's function (which is the ground motion for a unit constant slip in the subfaut k-th for the frequency f). Sk is the source function defined in the frequency domain whose temporal duration is given by the rise time Rk. Each subfault is allowed to slip once. The source function Sk is modulated by the slip amplitude Uk and is shifted in time by the rupture time tk. Further details on the method can be found in Hernandez et al. (1999). Because we calculate the complete Green's function we simulate the static displacement at zero frequency. This static ground motion Di, computed at each benchmark i, is used to model the ground deformation field resulting from DIn-vector d are related by the model operator g as d = g(m). The elements of d are real and imaginary parts of the displacement spectrum of the 3 components for selected stations band-pass filtered between 0.1 and 1.5 Hz. The g operator is a non-linear function of tk and Rk. If we assume an initial parameter vector m 0, we can get the solution mi+1 by an iterative approach. Using the iterative inversion algorithm based on the work of Tarantola (1987), the inverse solution, mi+1, is given by where b is a damping constant with values between 0 and 1, which is used to prevent divergence (we have chosen b = 0.05). A i is the Jacobian matrix of g(mi). All the derivatives are computed analytically in the frequency domain.

Fault parameterization
The goal of this study is to model geodetic and strong motion data collected during the three largest shocks of the 1997 Colfiorito seismic sequence. We first invert these data to image the rupture history of the two earthquakes of September 26th occurred at 00:33 and 09:40 UTC, respectively. The hypocenters of these two events are taken from Amato et al. (1998) and Deschamps et al. (2000): the first one (00:33 event) is located at the west corner (43.0305°N, 12.8622°E) at the bottom of the fault, the second (09:40) event (43.0255°N, 12.8917°E) is placed at the east fault corner at the bottom of the fault (as shown in fig. 4). The fault geometry is constrained from the CMT solutions proposed by Ekström et al. (1998): the dip and the strike of the 00:33 event are 46°and 152°, respectively, while for the 09:40 event they are 42°and 144°. The slip direction is consistent with a pure normal faulting with a rake angle equal to 270°. The fault that ruptured during the 00:33 event has both length and width equal to 7.5 km, while the 09:40 fault is 12.5 km long and 7.5 km wide. This latter fault is shallower than the former: the depth of the top of the fault is 0.7 km for the 09:40 rupture and 1.5 km for the 00:33 event. We parameterize the fault planes assuming 9 subfaults for the 00:33 fault and 15 subfaults for the 09:40 rupture plane. This parameterization is rather raw, but it is necessary to reduce the number of unknowns so that the problem does not become underdetermined. We simulate ground motion waveforms up to 1.5 Hz adopting a source time function for each subfault: we use a smooth ramp function described by three parameters: the final slip amplitude (U k), the rise time (Rk, corresponding to the local slip duration) and the rupture time (tk).
The fault geometry for the October 14th Sellano earthquake is taken from the model (Mod4) proposed by Salvi et al. (2000). The strike, dip and rake of the fault are equal to 135°, 45°and -90°respectively. The fault is 9 km long and 6 wide and it is divided into 24 square subfaults. The hypocenter (42.9190°N, 12.9260°E) is located 3 km along the strike from the northern edge of the fault and 4.5 km along the dip from the top of the fault, which is located at a depth of 2.4 km.

The two September 26th main shocks
In this study we invert the GPS measurements and the DInSAR observations in order to image the slip distribution for the two main shocks of September 26th (00:33 and 09:40 UTC, respectively). In this first inversion attempt we also slightly change the source geometry with respect to that described in the previous section to investigate how the fault position modifies the inverted source model. We compare the slip distribution resulting from our inverse approach with those proposed in the literature. We perform our inversions either by inverting each data set independently or by a joint inversion of GPS and DInSAR observations. The results of these calculations are shown in fig. 5a,b, while in table IV we list the fault parameters used in each source model that are modified with respect to those described in the previous section.
We associate an error estimate, introduced through the diagonal elements of matrix Cd, with each point of the 194 range change values extracted from DInSAR data. We have chosen a standard deviation equal to 3 cm for each DIn-SAR point, while for the GPS data we rely on the error estimates proposed by Hunstad et al. (1999). We assume that all the data are independent (off diagonal elements equal 0 in the Cd matrix). We also introduce a relative weight to all the data. The adopted weight values are a factor 1 for DInSAR points, a factor 10 for the horizontal components of 2 GPS stations for which the displacement is significantly higher than the error ellipse (stations CROC and PENN in fig. 3), a factor 5 for horizontal GPS points in the fault vicinity for which the displacement is slightly larger than the error ellipse (such as CAPA and FOLI) and a negligible weight (10 -5 ) for the GPS vertical component and the others GPS stations. This weighting is chosen in order to represent appropriately the intrinsic error of each data set. We first investigate the effect of changing the dip angle and the fault position with respect to the source model described above. We translate the fault in the direction perpendicular to the fault strike and we compare the misfit to the observed data to test the inversion results. We perform several least-square inversions using the fault parameterization described above and a homogeneous initial model m0 (slip is equal to 32 cm and 57 cm for the 00:33 and the 09:40 events, respectively). The model parameters and the retrieved final seismic moments for the different initial source models are listed in table IV: models M13 to M18 are suitable to test the influence of the assumed dip angle, while the others test the position of the two faults. The retrieved source models are plotted in fig. 5a,b for the 00:33 and the 09:40 earthquakes, respectively. For each initial model, we show the inverted models obtained either through the combined inversion of GPS and DInSAR or by the separate inversion of each distinct data set. We check the fit to the observed data by measuring the misfit (i.e., the difference between observed and syn- Table IV. Fault geometry and position for the source models shown in fig. 5a thetic waveforms) with the standard deviation and we list the score of each model and the available resolution in table IV. The best fits are obtained with models M7, M8 and M9; in this case the fault planes are shifted by 1.5 km toward the SW perpendicularly to the strike direction. All the inverted models shown in fig. 5a,b point to a common pattern: the two main shocks are composed of a main slip concentration located close to the nucleation point. For the 00:33 event the slip is lo-cated in the northern part of the fault near the hypocenter location and all the inversions show a single patch generally located at a depth between 4 and 6 km and no more than 6 km SE of the hypocenter. For the 09:40 event, the moment release is mainly located in the first 7 km NW of the focus. A second smaller slip concentration appears in the upper part of the fault at about 10 km North of the nucleation point in some inversion results. The inversion of GPS data alone has a smaller resolution than that of DInSAR observations, as expected due to the paucity of GPS benchmarks. The source models inferred by inverting GPS data require having slip at shallower depth than DInSAR observations. Dip angles between 37°and 47°do not substantially affect the inversion results. It is important to note that with the model geometry M1 the PENN GPS station is located on the hanging wall side of the fault while the measurements indicate a motion of this station in the northeast direction. This means that the fault is either located more south or deeper. The fit to the GPS displacements resulting from model M7 (which is obtained from the combined inversion of GPS and DInSAR data) is shown in fig. 3 while the fit to DInSAR observations is displayed in fig. 6a. This source model yields the best fit to the data. GPS displacements vectors are well reproduced although we underestimate the horizontal displacement at PENN. The vertical deformation resulting from GPS measurements is not used in the inversion, however model M7 provides a good fit to vertical displacement at the COLF and CROC benchmarks. We show in fig. 7 the comparison between observed GPS displacements and those simulated by using model M9 obtained by inverting GPS data only and has the same fault parameters of M7. The comparison between figs. 3 and 7 shows that model M9 reproduces better the horizontal displacements at PENN and CROC, but it produces a different orientation at CAPA. On the contrary, model M7 better fits the orientations of observed displacement vectors underestimating the amplitudes at PENN and CROC. These results confirm the finding of Hunstad et al. (1999), who concluded that slip at shallow depth is needed to fit the amplitude of the GPS displacement at PENN. Figure 6b shows the fit to DInSAR observations obtained through model M8. The comparison between fig. 6a,b shows that model M8 produces slightly smaller residuals than M7; in any case, these two models yield a similar fit to the observed data. Figure 8 shows the slip distribution for the two events of September 26th (model M7) and summarizes the fit to GPS and DInSAR by means of a histogram of the difference between data and simulated displacements. We will use this model, obtained through a combined inversion of GPS and DInSAR observations, to fit the strong motion data. However before performing any modeling of seismograms, we explored the space of acceptable solutions using a genetic algorithm. Using the fault geometry of model M7, the slip amplitude on each subfault is left free to vary between 0 and 1.5 m. The set of solutions found by the genetic algorithm is shown in fig. 9. The main goal of this genetic inversion is to obtain a set of acceptable models (with similar misfits) to have in-  Table 4. The two histograms to the bottom show count as a function of the residual difference between data and synthetics for each data set. For the GPS histogram the dark bars correspond to the GPS point having a significant weight in the inversion (horizontal components of benchmarks CROC, PENN, COLF CAPA and FOLI) and the light ones to the GPS data not used in the inversion.
formation about the uniqueness of the solution. It is interesting to note that, even if the misfit values obtained with the source models plotted in fig. 9 are similar, the genetic inversion yields a wide variety of acceptable models. The comparison among the models obtained through the genetic algorithm gives direct information on the available resolution of the slip distribution. For this reason we have plotted in fig. 9 the results of the genetic algorithm classified for three range of down-dip distance along the fault plane. This test emphasizes that the shallow portion of the fault plane is well resolved, with the exception of the north-western termination of the 09:40 fault (where the second slip concentration is located, as shown in fig. 5a,b and of the south-eastward termination of the 00:33 fault. On the contrary, the deepest portion of the two fault planes is not well resolved. This implies that it might be difficult to constrain with geodetic data the absolute values of slip and the details of slip gradient in the deepest portion of the fault. Figure 9 also shows that the result of the least square inversion (M7 is represented by the dashed lines) is in good agreement with the genetic algorithm and can be considered a reliable average among all the possible solutions. The maximum slip zone found with the genetic algorithm is indeed observed in the middle and the bottom of the fault close to the hypocenter location (SE part of the fault). For the 00:33 event the mirror pattern is observed. The maximum slip is located at the bottom on the western edge of the fault and the slip amplitude is better resolved in the northern part of the fault probably due to the poor data coverage in the south. Fig. 9. Slip profiles along the strike direction for three different depth ranges resulting from a genetic algorithm inversion of DInSAR and GPS data. 105 models were tested from a model space with about 1020 significantly different models using a population size of 200 and 500 generations. The best 20 models are plotted in thin black lines, the next 180 acceptable models found with the genetic algorithm are plotted in thin gray lines, the solution found with the leastsquares formulation is plotted in thick dashed lines. Fig. 10. Slip distribution, resolution of model parameters and fit to the data (top, middle and bottom panels, respectively) obtained from the least-squares inversion of the DInSAR data for the October, 14th Sellano earthquake. The fault is parameterized with 24 square subfaults (drawn on the fault plane), each of which has length of 1.5 km. The contour line interval for the slip distribution is 0.2 m, while for the resolution is 10%. The histograms show counts as a function of the residual difference between DInSAR data points and synthetics. data set on average. We have also verified the influence of the initial fault size and resolution by allowing the slip to vary over a larger fault area, by using 48 subfaults instead of 24. However, the final fit to the data is not improved, although the number of parameter is doubled.

Modeling strong motion waveforms through a forward and inverse approach
7.1. Modeling the strong motion waveforms recorded during the 00:33 event We use the slip distribution imaged for the 00:33 event by inverting DInSAR and GPS observations (model M7) to fit the strong motion waveforms through a forward approach. We assume a constant rupture velocity (2.6 km/s) and a constant rise time equal to 1 second: the resulting rupture model is shown in Fig. 11. The fit to the displacement time histories obtained from the accelerograms recorded at Assisi (ASS010), Matelica (MTL) and Cerreto di Spoleto (CERT) are shown in fig. 12. This figure points out that the slip model inferred by inverting geodetic data yields an acceptable fit to the observed ground motion waveforms. This source model underestimates the amplitudes observed on the horizontal components recorded at the

The October 14th Sellano event
The sparse configuration of the GPS network in the Sellano area (see figs. 1 and 3) did not allow us to measure significant static displacements (Salvi et al., 2000). Thus, the only available geodetic data to study the Sellano earthquake are the DInSAR observations. This data set is composed of 300 points extracted from the fringe pattern obtained from August 9th and October 17th images. It is not necessary to consider a larger data set because in this case the information redundancy would involve singularity problems using the least-square inversion scheme. Also in this case, we have chosen a standard deviation equal to 3 cm for each point. This estimated value is mainly ascribed to the ambiguous relationship between the phase difference measured with the satellite and the modeled absolute displacement.
The least-square inversion using this parameterization and a homogeneous initial model m 0 (slip is equal to 30 cm to have an initial moment consistent with the one found by Ekström et al., 1998) is shown in fig. 10. The slip is located in the north bottom part of the fault. The slip resolution is better at shallow depth probably due to the higher surface sensitivity to shallow slip (Hernandez et al., 1999). This model fits extremely well the DInSAR data. Less than one third of a fringe remains unexplained for all the Fig. 11. Rupture model used to simulate synthetic ground displacements for the September 26th, 00:33 UTC, earthquake. The slip distribution is that obtained by inverting geodetic data (model M7 of fig. 5a,b); the rupture time distribution is derived by assuming a constant rupture velocity equal to 2.6 km/s and rupture initiation at the hypocenter. The rise time is constant and equal to 1 s; contour lines are spaced 0.2 m and 1 s for slip and rupture times, respectively. mograms recorded during the main shock of the sequence (see fig. 14a). The fit to the other available stations is not so good (see fig. 14b), although the vertical component is sufficiently reproduced in amplitude with a modest misfit in phase.
We therefore invert the observed waveforms using the rupture model shown in fig. 13 as a starting model in the numerical approach previously described. We invert the observed ground displacements obtained from the recordings at the four stations shown in fig. 14a to constrain the slip and rupture times distribution on the fault plane. After several tests to find the best convergence, we chose a smaller a priori variance for the slip amplitude (Cm = 49) than for the rupture times (C m = 100). We assigned a greater weight to the digital stations AS010 and CERT, which are not contaminated by site effects and have an acquisition system with larger dynamic range, by giving a smaller variance (Cd = 1/16), while we use Cd = 0.25 for MTL and NCR. This latter station is particularly important since it is located on the hanging wall of the fault. The final slip distribution inferred for the 09:40 event is shown in fig. 15a. The main slip zone is still located on the southern portion of the fault but it is slightly shallower than that inferred inverting geodetic data (compare figs. 13 and 15a). A second asperity is 1370 Bruno Hernandez, Massimo Cocco, Fabrice Cotton, Salvatore Stramondo, Oona Scotti, Françoise Courboulex and Michel Campillo other stations (BVG, CLF NCR, CSA). This might be due to local site effects affecting the observed recordings; this is certainly true for Colfiorito (CLF) and Nocera Umbra (NCR). In this study we will not invert the ground motion waveforms to constrain the rupture history. Despite several attempts, we did not improve the fit to the recorded data obtained through the forward approach described above.

Waveform inversion for the 09:40 main shock
We start the study of the rupture process for the 09:40 earthquake performing a forward modeling of ground motion waveforms using the slip model which provides the best fit to both GPS and DInSAR data (model M7). We again assume a constant rupture velocity (2.6 km/s) and the nucleation at the southern edge on the fault bottom; the rise time is also constant and equal to 1 s. This rupture model is shown in fig. 13. We use this source model to compute synthetic ground displacements at selected sites: the fit to the recorded waveforms is shown in fig. 14. This test demonstrates that the slip model constrained by inverting geodetic data is able to fit quite well some of the seis-   13. Rupture model used to simulate synthetic ground displacements for the September 26th 09:40 UTC earthquake. The slip distribution is that obtained by inverting geodetic data (model M7 of fig. 5a,b); the rupture time distribution is derived by assuming a constant rupture velocity equal to 2.6 km/s and rupture initiation at the hypocenter. The rise time is constant and equal to 1 s; contour lines are spaced 0.2 m and 1 s for slip and rupture times, respectively.  fig. 13. In a) we show the fit to the observed time histories that will be used in the waveform inversion; b) shows the fit to those waveforms that will be not used. a b still present in the shallow northern section of the fault. The rupture time distribution is different from that assumed for the starting model and reveals high rupture velocity at shallow depths. The fit to the observed data is shown in fig. 15b for the four stations used in the inversion. The total variance reduction for this inversion is 37% and the solution is found after 21 iterations. The fit to the other recorded data is not improved with respect to that shown in fig. 14b and obtained through the forward approach. The final seismic moment is 8.2 × 10 17 Nm, which is slightly less than the values inferred by modeling regional and teleseismic data (M0 = 12.0 × 10 17 Nm: Ekström et al., 1998 andM0 = 11.4 × 10 17 : Dziewonsky et al., 1998).

Modeling strong motion waveforms recorded during the October 14th Sellano earthquake
Also for the October 14th Sellano earthquake we start our analysis performing a forward modeling of recorded strong motion waveforms using the geodetic slip model and assuming a constant rupture velocity equal to 2.6 km/s and rise time equal to 1 s. The resulting source model is shown in fig. 16a. The fit to recorded ground displacements is shown in fig. 16b,c. This figure shows that the fit to recorded data is not very good, because this source model overestimates the ground motion amplitudes recorded at most of the stations (such as Assisi -AS010 -and Foligno -FOL). We therefore conclude that for the Sellano event, the slip distribution inferred from the inversion of geodetic data and a constant rupture velocity do not fit the recorded seismograms.
Thus, we performed a waveform inversion attempt: we inverted the three components of the ground displacements recorded at ANN, COL, AS010, and CERT. All these stations were equipped with digital acquisition systems and have recorded the entire radiated waveform (P and S waves). In our inversion we gave a larger weight to the station CERT located south of the event (Cd = 1/16.0) and less weight to the others stations located north of the event (Cd = 1/9.0). The inversion results are displayed on fig. 17a and the fit in time between Fig. 15a,b. a) Slip and rupture time distributions resulting from the iterative inversion of ground displacements for the September 26th 09:40 UTC event. The rupture model shown in fig. 13 was used as starting model in the iterative approach. The rise time is assumed to be constant (1 s). b) Comparison between recorded (solid line) and synthetic (dashed line) ground displacements for the 09:40 earthquake calculated using the source model resulting from the waveform inversion and displayed in (a).   fig. 17b,c. The total variance reduction for this inversion is 49% and this solution, which matches the data much better than the geodetic model, is found after 35 iterations.

Summary of the results and conclusive remarks
Despite their moderate magnitudes, the three largest shocks of the 1997 Umbria-Marche seismic sequence allowed the collection of a quite complete data set. Coseismic ground deformations were measured by GPS displacements and DInSAR interferograms, while strong motion data were recorded near the causative faults. The analysis of the aftershock sequence (Chiaraluce et al., 2003; and the analysis of broad band regional data (Morelli et al., 2000;Pino and Mazza, 2000) provided reliable information on the earthquake locations, the focal mechanisms and the fault geometry. In this study we have used these data and the results of previous studies to constrain the fault geometry and to image the rupture history for the three main shocks that occurred on 26th September 1997 [at 00:33 and 09:40 UTC, respectively] near Colfiorito and on 14th October 1997 near Sellano. The modeling of both geodetic measurements and near-field strong ground motion data, either through a forward or an inverse approach, yielded the following results: i) The 26th September 1997 earthquake that struck the Colfiorito area at 00:33 UTC is composed of a single slip patch located close to the hypocenter at depth ranging between 3 and 6 km; the rupture propagation is unilateral toward the south-east direction. It ruptured a normal fault striking along the Apennines and dipping 46°to the SW.
ii) The second 26th September 1997 event that struck again the Colfiorito area at 09:40 UTC has broken in the opposite direction (toward the northwest) a distinct normal fault dipping 42°to the SW. The inferred source model for this event consists of a larger slip patch located close to the hypocenter and a second smaller slip concentration, shallower than the other, and located at the northeastern edge of the fault (nearly 10 km from the hypocenter). This model is able to fit both geodetic data and ground motion waveforms.
iii) The 14th October 1997 Sellano earthquake also consists of a single slip patch located at the bottom of the fault in its northeastern portion. iv) These events are characterized by evident rupture directivity effects. Our results are in agreement with previous investigations (Capuano et al., 2000;Pino and Mazza, 2000). Observed ground displacements can be reasonably fitted with a constant rupture velocity model (we use in this study 2.6 km/s). However, the inverted rupture models show a more heterogeneous rupture time distributions, which are consistent with the directivity effects. v) For all the three main shocks the slip is concentrated at depth ranging between 4 and 6 km. The lack of coseismic slip at top of the fault (at shallow depth) suggests that the observed surface breaks are caused by induced deformation, consistently with some interpretations of field observations (Cinti et al., 1999). We emphasize here that the source models obtained by inverting geodetic data have an excellent resolution of the shallower portion of the fault planes.
The main goal of this paper is to constrain the slip distribution during the main shocks of the Umbria-Marche sequence. Although many recorded seismograms are affected by local site amplification effects, the waveform modeling performed in this study, either through a forward or an inverse approach, yielded an acceptable fit to recorded waveforms. We believe that the quality of the match between synthetic and observed seismograms is quite good as can reasonably be expected, considering the uncertainties in the crustal velocity model and given the lack of absolute time of the recordings. It is important to point out that the inferred source models allow the fit to both geodetic (GPS and DIn-SAR) and ground motion waveforms.
In our opinion, there are two further issues that are important to discuss: the first concerns the position of the faults resulting from our best fitting model, the second involves the possibility to constrain the rupture velocity. The best fit to available observations collected for the two main shocks of September 26th is obtained for two fault planes whose position is shifted by 1.5 km to the SW perpendicularly to the strike direction inferred from the CMT focal mechanisms. The size of this spatial translation lies within the location error estimated for these two seismic events (Amato et al., 1998;Chiaraluce et al., 2004). Thus, we believe that the positions of the fault planes inferred in this study cannot be considered very different from those proposed by Chiaraluce et al. (2003).
The rupture time distributions imaged through our waveform inversion procedure are very heterogeneous and characterized by a variable rupture velocity. However, we were able to reproduce sufficiently well the observed recordings also by means of a source model characterized by a constant rupture velocity. Thus, our modeling results demonstrate that we can model the observed time histories of ground displacements with very different rupture history models. Some of the inferred models show super-shear rupture velocity in specific fault patches. We cannot exclude this possibility. However, due to the low frequency bandwidth investigated in this study, the small dimensions of the fault planes and the non-uniqueness of the inverse approach we cannot provide a well constrained rupture history model. Further investigations are needed to unravel this problem.

Introduction
In the estimate of seismic risk, the determination of ground motion parameters like spectral characteristics, peak ground displacement and peak ground acceleration is very important for a quantitative assessment of the problem. The knowledge of these parameters provides the basis for a classification of the territory and identifies the areas for which great damage in case of strong earthquakes is expected. The determination of ground motion parameters can be carried on by numerically simulating the time history related to the maximum expected earthquake in a given area. Among the techniques used for the estimate of ground motion parameters, the stochastic method (Boore, 2003) is widely used to predict the ground motion due to a seismic input, which can be modeled taking into account source, path and site effects. The method provides the ground motion by combining a parametric description of amplitude spectrum with a random phase spectrum. This technique has been applied to estimate peak parameters of ground motion in many regions of the world (Akinci et al., 2001;Castro et al., 2001Atkinson et. al., 2002 and is particularly useful to simulate the ground motion for the frequency range usually investigated by earthquake engineering.