Vol51,2_3,2008 461 ANNALS OF GEOPHYSICS, VOL. 51, N. 2/3, April/June 2008 Key words data inversion – geodesy – coseismic dis- placement 1. Introduction The 1997 Umbria-Marche sequence in- cludes two moderate magnitude earthquakes (September 26th, Mw = 5.8 and Mw = 6.0) that occurred near the village of Colfiorito followed by a sequence of aftershocks lasted several months. The aftershock activity included a Mw = 5.6 earthquake (that occurred near the village of Sellano on October 14th) and three events with magnitude larger than 5.2. The Colfiorito sequence has been extensively studied using ge- ologic data, strong motion data, teleseismic and regional waveforms, GPS, levelling and SAR data (see e.g. Cinti et al., 1999; Capuano et al., 2000; Pino et al., 1999; Hunstad et al., 1999; Lundgren and Stramondo, 2002; De Martini et al., 2003; Crippa et al., 2006). Larger shocks share SW strikes and normal fault mechanisms. After the September 26th event, 12 monu- ments from the IGMI95 GPS network within 30 km off the epicentre were reoccupied, including FOLI and COLF that were occupied continu- ously during the survey. Coseismic diplace- ments were obtained by comparing the 1997 survey coordinates with those collected in 1995. Two sites (CROC, -13.8 cm N and -1.9 cm E; PENN, 7.6 cm N and 8.0 cm E) showed Inversion of synthetic geodetic data for the 1997 Colfiorito events: clues on the effects of layering, assessment of model parameter PDFs, and model selection criteria Antonella Amoruso and Luca Crescentini Dipartimento di Fisica «E.R. Caianiello», Università degli Studi di Salerno, Italy Abstract The 1997 September-October Umbria-Marche sequence has been extensively studied in the past by analyzing coseismic displacement data (GPS, leveling, SAR). Here we focus on synthetic data representative of the main event of the 1997 Umbria-Marche sequence and investigate the effects of a crustal layering proper to the Colfior- ito area on surface displacements and inferred source features when inverting coseismic geodetic data without taking into account layering. We compare bootstrapping and NA-Bayes as tools for parameter uncertainty as- sessment and show how the Akaike Information Criterion can be used to select the model which is most likely to be correct. Since SAR images offer the most complete coverage of the study area, we use synthetic line-of- sight displacement data. Mailing address: Dr. Antonella Amoruso, Dipartimen- to di Fisica «E.R. Caianiello», Università degli Studi di Sa- lerno, Via S. Allende, 84081 Baronissi (SA), Italy; e-mail: antonella.amoruso@sa.infn.it Vol51,2_3,2008 4-03-2009 10:28 Pagina 461 462 A. Amoruso and L. Crescentini significant horizontal displacements and one (CROC, -24.7 cm) showed significant vertical displacements (Hunstad et al., 1999). The epicentral area is crossed by Line 21 of the Italian first-order levelling network (operat- ed by IGM), which runs from Foligno, central Umbria, to Fiumesino, on the Marche coast; the first (southernmost) section of Line 21 trends almost N-S between Foligno and Fossato di Vi- co and includes about 40 bench marks. The route was measured in 1951, 1992 and was re- surveyed in 1998. The data showed maximum subsidence of 7.8 cm, and maximum uplift of 1.6 cm (De Martini et al., 2003). Differential interferograms related to the 1997 Colfiorito events have been obtained from ERS-1 and ERS-2. Because of the strong topo- graphic relief of the Apennines and the signifi- cant vegetation cover of the epicentral area, even the best interferograms show large decor- relation areas (Salvi et al., 2000). For this rea- son, even though a large number of SAR data are available from July 1993 to October 1997, only few significant interferograms have been obtained, covering different periods of the se- quence. The most studied interferogram is the 35-day ERS-2 one (e.g. Stramondo et al., 1999; Crippa et al., 2006) covering the period 7 Sep- tember to 12 October 1997. Geodetic data are usually affected by both correlated and uncorrelated noise. For example, levelling data are obtained by measuring the el- evation difference between consecutive bench marks (section height differences). Bench mark elevations with respect to the reference bench mark can then be computed, and are often used instead of section height differences for model- ing. This operation introduces nondiagonal terms in the covariance matrix of bench mark elevation differences, whose diagonal terms are related to uncorrelated noise (e.g. Amoruso and Crescentini, 2007). As regards SAR data, the unwrapped differ- ential phase is also affected by both uncorrelat- ed and correlated noise. Correlated noise is caused by atmospheric disturbances, and can be expressed by an exponential autocorrelation function. This noise can be handled following two different schemes: using a non-diagonal covariance matrix in the data inversion proce- dure (Fukushima et al., 2005) or removing from the interferogram the influence of atmospheric heterogeneities after considering stable areas close to the epicentral area to analyse the spa- tial autocorrelation (e.g. Crosetto et al., 2002). We prefer the former scheme, since, as stressed in Crosetto et al. (2002) the removal procedure «suffers for its intrinsic limitation related to the non-stationarity of the analysed signal, which practically confines its use to the vicinity of the stable areas». To give an idea of overall noise in SAR data related to the Colfiorito events, Lund- gren and Stramondo (2002) assigned a standard deviation of 1 cm to each contoured data point (i.e., points of line-of-sight, LOS, displacement contours at 28 mm intervals). If both correlated and uncorrelated noise is present, measured data are not statistically in- dependent. However, the data covariance ma- trix is symmetric and can be reduced to diago- nal form by means of a rotation matrix; the same rotation transforms the data to independ- ent form (i.e. a set of statistically independent data, which are linear combinations of the measured ones). The eigenvalues of the covari- ance matrix give the uncertainties of the trans- formed independent data. This procedure can also be used in case of non-normal distribution of errors and consequently for robust fitting (Amoruso and Crescentini, 2007). Model parameters are usually obtained from experimental data through the appraisal of a cost function which measures the disagreement between data and model predictions. Several different causes contribute to differ- ences between predictions and observations (the residuals of the retrieved model): measure- ment errors, intrinsic misfits (e.g. misfits due to heterogeneities in slip distribution, non-planar faults, heterogeneities in the elastic properties of the medium, when the adopted model is too simple, like a planar uniform-slipping fault in a homogeneous half-space), local response ef- fects (e.g. soil compaction, landslides, topo- graphic effects). For example, Amoruso et al. (2004) generated synthetic geodetic coseismic data (GPS, SAR, and levellings) in a layered medium and performed inversions of the syn- thetics, under the assumption of homogeneous half-space. Horizontal displacements are more Vol51,2_3,2008 4-03-2009 10:28 Pagina 462 463 Inversion of synthetic geodetic data for the 1997 Colfiorito events: clues on the effects of layering, assessment of model parameter PDFs, and model selection criteria affected than vertical ones by the presence of a superficial soft layer, but differences with syn- thetics generated in a homogeneous half-space did not show any simple regular behaviour, even if a few features could be identified. Con- sequently, also retrieved parameters of the ho- mogeneous equivalent fault obtained by uncon- strained inversion of surface displacements did not show a simple regular behaviour. Amoruso et al. (2004) pointed out that the presence of a superficial layer may lead to misestimating sev- eral fault parameters both using joint and sepa- rate inversions of the three components of syn- thetic displacement. The effects of the presence of the superficial layer depend on whether all fault parameters are left free in the inversions or some of them are fixed a priori. In the inversion of any kind of coseismic geodetic data, fault size and slip could be largely misestimated, but the seismic potency (product of ruptured area and average slip) is well determined (within a few per cent) in most cases even neglecting lay- ering. Of course, the amount of the effect of the superficial low-rigidity layer depend on the contrast of the elastic parameters. For example, Battaglia and Segall (2004) and Crescentini and Amoruso (2007) have shown that layering ef- fects are negligible and relevant respectively when considering two different volcanic envi- ronments and layering properties (Long Valley and Campi Flegrei). Assessing the effects of layering is a case-to-case problem. Intrinsic misfits and local response effects, as well as possible systematic errors, could be spatially coherent over small or large areas, and are not zero-mean random variables. Rigorous handling of data should take into account, and properly treat, all error sources (measurement errors, systematic errors, intrinsic misfits and so on), after a reliable a priori estimate of all the errors, and the determination of how much dif- ferent error sources effectively affect data. This is quite a hard task, largely case-to-case de- pendent. In case of levelling data, if measure- ment errors can be reliably estimated, after in- verting the transformed independent data, re- duced chi-square of residuals can be computed and used to refine uncorrelated error variance when the number of degrees of freedom is eas- ily computable, e.g. in case of one uniform-slip rectangular fault in a homogeneous half-space (Amoruso and Crescentini, 2007). The cost function is usually obtained from maximum-likelihood arguments according to the assumed statistical distribution of the resid- uals. Misfit is written in terms of the trans- formed independent data set, and, usually, ei- ther the mean squared deviation of residuals (chi-square fitting, M2, proper for normally dis- tributed errors) or the mean absolute deviation of residuals M1, proper for two-sided-exponen- tially distributed errors and robust fitting) are used. If possible, it is preferable to compare re- sults obtained using M1 and M2. Estimates of the model parameters can be ob- tained by the best-fitting approach, i.e. the mini- mization of the cost function. When inverting ge- odetic data, the cost function depends nonlinear- ly on parameters and is characterized by a rough landscape and several local minima. Monte Car- lo techniques are the only efficient tool to search for the global minimum, but require computing a large number of forward models (data prediction given a set of model parameters). Thus the use of fast codes for forward mod- elling is preferable if not essential. Among the Monte Carlo global optimization methods (e.g. uniform sampling, genetic algorithms, simulat- ed annealing) we prefer to use Adaptive Simu- lating Annealing (ASA, Ingber 1993). General- ly speaking, each step of any Simulated-An- nealing algorithm replaces the current point in the parameter space by a random «nearby» point, chosen with a probability that depends on the difference between the corresponding cost- function values and on a global parameter (T) that is gradually decreased during the process according to some annealing schedule. ASA considers that different parameters might re- quire different annealing schedules. The expo- nential annealing schedules used by ASA en- sure ample global searching in the first phases of search and ample quick convergence in the final phases. Even if success in finding the global minimum is never guaranteed and inver- sions from different starting points and/or using different program options can lead to different endpoints, the ability of ASA to self optimize its program options recursively makes it a very powerful technique. Vol51,2_3,2008 4-03-2009 10:28 Pagina 463 464 A. Amoruso and L. Crescentini The best-fitting approach does not give pa- rameter uncertainties and consequently is of limited significance. A commonly used tech- nique (associated with global minimization) to assess acceptable parameter ranges consists in the use of the bootstrap percentile method. This method applies the best-fit technique to a large number of synthetic data sets, obtained by ran- domly resampling from the actual data set with replacement and is often used to estimate confi- dence intervals without making assumptions about the underlying statistics of errors. Bootstrapping also estimates correlations between different parameters, by forming a scat- ter plot of all the estimates for each parameter pair and visualizing the correlations between the parameter pairs (e.g. Amoruso et al., 2005). The computed values for the model parameters from the synthetic data sets form an estimate of the sampling distribution of the parameter values, independently from the underlying statistics. This method gives reliable results only if the hy- pothesis of independent identically-distributed data is verified. As a consequence, in the pres- ence of a non-diagonal data covariance matrix, bootstrapping can only be applied to the inde- pendent rotated data, not the measured data. A different approach is used in the Neigh- bourhood Algorithm (NA-sampler, Sambridge, 1999a). NA generates ensembles of models which preferentially sample the good data-fit- ting regions of the parameter space, rather than seeking a single optimal model. The algorithm makes use of only the rank of a data fit criteri- on rather than the numerical value; in other words, points in the parameter space are listed according to their capability to fit the data (an- swering the question «does point A give a bet- ter fit to the data than point B?») without quan- tifying the difference in a precise way. With this rank-based approach, the weight of each of the previous sampled points in driving the search depends only on their position in the rank list. The companion code NA-Bayes (NAB, Sam- bridge 1999b) consists of an algorithm for us- ing the entire ensemble of models produced by NA, and deriving information from them in the form of Bayesian measures of resolution, co- variance and marginal probability density func- tions (PDF) etc. The input ensemble in princi- ple can follow any distribution and be generat- ed by any search method (e.g. the NA-sampler algorithm, a genetic algorithm or simulated an- nealing algorithm). In any case, we assume a priori a physical model to explain the data. In principle, any set of data can be almost perfectly fit by using a sufficiently complicated model, but it could be unrealistic and overmodelling (i.e., using more parameters than necessary) has to be avoided. Attempts to find the model that best explains the data with a minimum number of parameters accomplish the parsimony principle. If nested models (the more complicated model includes the simpler one as a particular case) have to be compared, the F-test approach could be used. The F-test is based on tradition- al hypothesis testing. Under the normal distri- bution hyphothesis, F-test compares the differ- ences in χ2 between nested models with the dif- ference expected by chance, i.e. it quantifies the relationship between the improvement in χ2 and the relative increase in degrees of freedom g. The null hypothesis is that the simpler mod- el is correct, i.e. that none of the additional pa- rameters is significantly correlated with the model. Unfortunately, the test is unreliable if residuals depart from the normal distribution even slightly. Moreover, the F-test is not valid if the models are not nested or it is difficult to es- timate the actual number of g (as in the case of distributed-slip faults, due to the presence of constraints like smoothness or slip positivity). The Akaike Information Criterion (AIC, Akaike, 1974) is a measure of the goodness of fit of an estimated statistical model. It is based on information theory and does not use hypoth- esis testing, so there is no conclusion about sta- tistical significance and rejection of a model. AIC attempts to find the model that best ex- plains the data with a minimum number of pa- rameters taking into account both the value of the likelihood function and the number of pa- rameters in the model, i.e. trading off the com- plexity of an estimated model against how well the model fits the data. AIC includes a penalty that is an increasing function of the number of estimated parameters. The model with the smallest AIC is most likely to be correct. For a small number of observations, a second order Vol51,2_3,2008 4-03-2009 10:28 Pagina 464 465 Inversion of synthetic geodetic data for the 1997 Colfiorito events: clues on the effects of layering, assessment of model parameter PDFs, and model selection criteria correction is added to the general AIC; correct- ed AIC (AICc) converges to AIC as the number of observations increases (e.g. Hurvich and Tsai, 1989). In addition, the difficulty of esti- mating the actual number of degrees of freedom can be surmounted by using the Akaike Infor- mation Criterion. Even if the results are not au- tomatically generalizable, Amoruso and Cres- centini (2007) have shown, in case of levelling measurements and the one-fault model of the 1908 Messina earthquake, that the Akaike In- formation Criterion is not only an effective tool to discriminate if increasing the number of sub- faults in a distributed-slip model really im- proves the model or does not, but also to esti- mate uncorrelated errors, which partially reflect model inadequacy. Here we investigate the effects of a crustal layering proper to the Colfiorito area on surface displacements and inferred source features when inverting synthetic geodetic data repre- sentative of the main event of the 1997 Umbria- Marche sequence. We compare different meth- ods for parameter uncertainty assessment and we show how the AICc can be used to select the distributed-slip model which is the most likely to be correct as a result of the inversion of avail- able data. Since SAR images offer the most complete coverage of the study area, we use synthetic LOS displacement data. 2. Generation and inversion of synthetic co- seismic displacement data 2.1. Generation To generate synthetic displacements we use the slip distribution in fig. 1, which resembles the recently published distributed-slip model of the main event in Crippa et al. (2006). The source fault (divided into 100-along-strike and 80-along-dip subfaults) is embedded both in a homogeneous half-space (HHS; Okada, 1985) 0 1 m 0 1 2 3 4 5 6 7 8 A lo n g − d ip d is ta n ce (k m ) 0 1 2 3 4 5 6 7 8 9 10 Along−strike distance (km) Fig. 1. Slip distribution on the source fault. Vol51,2_3,2008 4-03-2009 10:28 Pagina 465 466 A. Amoruso and L. Crescentini and the layered half-space (LHS; Wang et al., 2006) detailed in table I (Megna, 2007). Strike, dip and rake are 138°, 45°, and -75° respective- ly. We then add uncorrelated random noise. As previously mentioned (see Introduction) atmos- pheric disturbances generate correlated random noise that can be expressed as an exponential covariance function and would require trans- forming the data to independent form, but treat- ment of the effects of atmospheric disturbances is beyond the aims of this paper. Displacements (northward, eastward, vertical) were computed on a grid spanning 15000 m toward North, step 100 m, and 20000 m toward East, step 100 m. We simulated SAR-like unwrapped data through a proper combination of the displace- ment components taking into account the ERS- 2 LOS. LOS displacement (positive toward the satellite) maps for HHS and LHS are in fig. 2a and fig. 2b respectively. Figure 2c shows the residuals of the linear least squares regression between LHS and HHS displacements. As ex- pected (e.g. Amoruso et al., 2004) residuals are about 15% of the signal, and, being obtained using a linear best fit, are not merely due to a translation or a scale factor. LHS LOS displacement after the addition of noise (generated using random variates from a normal distribution with zero-mean and 1-cm- standard deviation) is shown in fig. 3. To reduce the number of data in the inversions, the dis- placement map was decimated (undersampling interval = 3 bins). Even if noise added to synthetics was ob- tained from a normal distribution, model inade- quacy is expected to enlarge tails of the residual distribution with respect to the normal distribu- tion and we minimize the mean absolute devia- tion of residuals M1, commonly used for robust fitting: (2.1) where n is the number of data and vi is the co- seismic change estimated for element i of the data set. The data set must be composed of sta- tistically independent data: if they are not (i.e. the covariance matrix is not diagonal) then they have to be transformed accordingly. Here ex- perimental data are not rotated because of the absence of correlated noise. Model predictions v–i (a) depend on model parameters a, and wi are the uncertainties of the independent data. 2.2. Inversion, uniform-slip model At first, we invert noisy LHS LOS displace- ments for a uniform-slip rectangular fault em- bedded both in HHS and in LHS using ASA. This model involves a parameter set a which consists of the strike and dip angles of the fault, the rake angle of the slip vector, two geometri- cal fault dimensions (along-strike length and along-dip width), location of start point of the fault upper side along the strike direction (x, positive Northward; y, positive Eastward), depth of the upper side of the fault (z, positive downward), and magnitude of the slip. Figure 4 shows the two best-fit uniform-slip faults over the distributed-slip source fault. Po- tency of the best-fit faults is the same inside 2% (� 2.65 × 107 m3), about 9% lower than the po- tency of the source fault (2.89 × 107 m3). The ( )a M w v v 1 1 i i i i n = - = r/ Table I. Multilayered model used in this work (Megna, 2007). Depth (km) VP (m/s) VS (m/s) Density (kg/m3) 0-1 2300 1200 2300 1-2 5400 2500 2600 2-8 6300 3400 2800 >8 5500 3000 2840 Vol51,2_3,2008 4-03-2009 10:28 Pagina 466 467 Inversion of synthetic geodetic data for the 1997 Colfiorito events: clues on the effects of layering, assessment of model parameter PDFs, and model selection criteria only noticeable difference between the two best-fit faults consists in the expected focusing effect in the HHS: while the LHS fault is cen- tered on the slip patch of the source fault, the HHS fault is shifted by about 700 m upward along the dip direction, i.e. 500 m in depth. PDFs of model parameters were estimated using both bootstrapping and NAB. Strictly speaking, a full optimization procedure should be performed at each resample in the bootstrap procedure, but in our case computation time would be unacceptable. In this work we use a faster method (downhill simplex; Press et al., 1992), starting from the ASA best-fit model to explore the parameter region surrounding it. As regards NAB, the input ensemble has been gen- erated by NA. Figure 5 shows estimated PDFs as well as best-fit parameter values. The true Fig. 2. Map of synthetic LOS displacements generated in the HHS (left) and LHS (middle). Right, map of the residuals of the best-fit linear regression between LOS displacements in the HHS and LHS. Fig. 3. Map of noisy synthetic LOS displacements generated in the LHS. Horizontal projection of the source fault in black. −29 5 cm 0 5 10 15 km 0 5 10 15 20 km Noisy LOS displacement −26 1 cm 0 5 10 15 km 0 5 10 15 20 km LOS displacement (homogeneous) −27 2 cm 0 5 10 15 km 0 5 10 15 20 km LOS displacement (layered) −15 29 mm 0 5 10 15 km 0 5 10 15 20 km LOS displacement (difference) Vol51,2_3,2008 4-03-2009 10:28 Pagina 467 468 A. Amoruso and L. Crescentini source values are also shown when pertinent, i.e. for strike, dip, and rake. The best-fit values and PDFs are in agreement with the true values, but strike is somewhat overestimated, owing to the uniform-slip approximation, as it will be clarified below. PDFs obtained using the two techniques are generally in mutual agreement; apparently larger differences relate to depth of the fault upper side and along-dip width in the LHS case. However, fig. 6 shows that there is a trade-off between them and the trade-off pattern is very consistent between bootstrapping and NAB. 2.3. Inversion, distributed-slip models As a second step, we consider a small num- ber of coplanar subfaults (2 along strike and 2 along dip, numbered along the strike direction, starting from the shallower part of the fault), which slip independently at the same rake an- gle, thus adding three additional parameters to the model. Figure 7 shows the two best-fit 2 × 2 faults over the distributed-slip source fault. Po- tency of the best-fit faults is the same inside 4% (� 2.6 × 107 m3), about 10% lower than the po- tency of the source fault (2.89 × 107 m3). As in the case of the uniform-slip model, the only no- ticeable difference between the two best-fit faults consists in the focusing effect. Once again, PDFs of model parameters were estimated using both bootstrapping and NAB, as detailed in the case of the uniform-slip fault. Figure 8 shows estimated PDFs as well as best-fit parameter values and suggests similar comments with respect to the uniform-slip model. PDFs from NAB are more scattered Fig. 4. Best fit uniform-slip fault in the HHS (thin-line rectangle) and LHS (tick-line rectangle) over the dis- tributed-slip source fault. Numbers give slip. 0.84 0.77 0 1 m 0 1 2 3 4 5 6 7 8 9 A lo n g − d ip d is ta n ce (k m ) 0 1 2 3 4 5 6 7 8 9 10 11 Along−strike distance (km) Vol51,2_3,2008 4-03-2009 10:28 Pagina 468 469 Inversion of synthetic geodetic data for the 1997 Colfiorito events: clues on the effects of layering, assessment of model parameter PDFs, and model selection criteria than from bootstrapping, especially for subfault slips. This behaviour is probably related to the cost function roughness: on the one side NA could be in some trouble in sampling the pa- rameter space efficiently, as suggested by its best-fit misfit (larger than that obtained by ASA, using the same number of forward mod- els), on the other, bootstrapping uses simplex, which could be unable to escape from local minima. This problem deserves great attention and deep investigation while inverting real data. As a first simple application of model selec- tion criteria (selecting the model that best ex- plain the data with a minimum number of pa- rameters) we use AIC corrected for small sam- ple sizes (Hurvich and Tsai, 1989) (2.2) to compare the uniform-slip fault model and the 2x2 one. Here n is the number of data, k is the number of parameters, and M is the chosen dis- crepancy (measure of lack of fit, e.g. Zucchini, 2000); in this case M = M1. Difference in AICc 1 2 lnAICc n k nk n n M= - - + 20500 21000 21500 x (m) 12000 13000 14000 y (m) 1500 2000 2500 z (m) 5000 5500 6000 6500 Along-strike length (km) 5500 6000 6500 Along-dip width (km) 138 140 142 144 146 Strike (deg) 43 44 45 46 47 Dip (deg) -80 -75 -70 Rake (deg) 0,7 0,8 0,9 1 Slip (m) Fig. 5. PDFs of inferred fault parameters in the HHS (thin lines) and LHS (thick lines): dashed lines, bootstrap- ping; solid lines, NAB. Vertical dotted lines show ASA best-fit values; vertical dash-dotted lines indicate true parameter values. Uniform-slip model. Vol51,2_3,2008 4-03-2009 10:28 Pagina 469 0 1 m 0 1 2 3 4 5 6 7 8 9 A lo n g − d ip d is ta n ce (k m ) 0 1 2 3 4 5 6 7 8 9 10 11 Along−strike distance (km) 0.23 0.59 0.46 1.43 0.17 0.57 0.56 1.44 470 A. Amoruso and L. Crescentini 5200 5600 6000 6400 6800 A lo n g − d ip w id th 1800 2000 2200 2400 z (m) 2D PDF Fig. 6. 2D PDF of fault upper-side depth and along-dip width as obtained using NAB (linear gray scale) and scatter plot from bootstrapping (circles). LHS uniform-slip model. Fig. 7. Best fit 2 × 2 fault in the HHS (thin-line rectangles) and LHS (thick-line rectangles) over the distrib- uted-slip source fault. Numbers give slip. Vol51,2_3,2008 4-03-2009 10:28 Pagina 470 471 Inversion of synthetic geodetic data for the 1997 Colfiorito events: clues on the effects of layering, assessment of model parameter PDFs, and model selection criteria is about -500 both for the HHS and the LHS, in- dicating that the uniform-slip model is too sim- ple and requires complex models to be tested. Increasing the number of subfaults would disclose the main features of the slip pattern and possible changes in fault geometry while relaxing the assumption of uniform slip (e.g. Amoruso et al., 2002), but also increases the difficulty of studying the cost function features. In this work, fault plane geometry and rake an- gle show very small changes when turning from the uniform-slip fault model to the 2 × 2 one, thus we use both the geometry and mechanism of the fault plane of the best fitting uniform-slip model to estimate the slip distribution, but we increase the along-strike length and the downdip width of the fault. To enlighten the ef- fects of neglecting layering when inverting ge- odetic data, we consider the HHS case. The fault is divided into an even grid of p-along- 21000 22000 23000 x (m) 10000 12000 14000 y (m) 1000 1500 2000 2500 3000 z (m) 4000 6000 8000 Along-strike length (km) 4000 5000 6000 Along-dip width (km) 130 135 140 145 150 Strike (deg) 44 46 48 50 52 Dip (deg) -85 -80 -75 -70 Rake (deg) 0 0,5 1 1,5 2 2,5 Slip #1 (m) 0 0,5 1 1,5 2 2,5 Slip #2 (m) 0 0,5 1 1,5 2 2,5 Slip #3 (m) 0 0,5 1 1,5 2 2,5 Slip #4 (m) Fig. 8. PDFs of inferred fault parameters in the HHS (thin lines) and LHS (thick lines): dashed lines, bootstrap- ping; solid lines, NAB. Vertical dotted lines show ASA best-fit values; vertical dash-dotted lines indicate true parameter values. 2 × 2 model. Vol51,2_3,2008 4-03-2009 10:28 Pagina 471 472 A. Amoruso and L. Crescentini strike x q-along-dip subfaults, with variable slip magnitude and constant rake. We minimize the cost function (2.3) using the algorithm for the least squares prob- lem with linear inequality constraints of Law- son and Hanson (1995). Here ∇2s is a finite dif- ference approximation to the Laplacian of the slip distribution and A is the area of each sub- fault. It would also be interesting to use misfit M1 of LOS displacements, but it leads to the prohibitive problem of nonlinear minimization in a huge space. Finding the optimal model parameters thus requires the simultaneous optimization of two objective functions, namely the fit to the data 1a w v v fA p q s 1 2 2 2 1 i i i i n i i p q # d - + # = = r ] c ] g m g/ / and the roughness, whose relative importance is controlled by the adimensional smoothing pa- rameter f. The choice of f is a basic problem in inverse theory. Traditionally, it is selected from a «trade-off» curve which plots the trade-off be- tween the two objective functions when varying the value of f. The trade-off curve is monotonic and asymptotes to the minimum value of the roughness (which is zero) at one end and the minimum value of the fit to the data at the oth- er. Choosing a suitable value of f is somewhat arbitrary, but the knee of the trade-off curve is a good compromise between the model complex- ity and the data fit. Cross validation is a more rigorous criterion (e.g., Matthews and Segall, 1993) but is computationally intensive. Tests performed by Arnadottir and Segall (1994) show that cross validation gives an optimal es- timate of f, while the trade-off curve gives a slightly smoother solution. We use the trade-off 1× 10 8 1× 10 9 1× 10 10 3× 10 10 7× 10 10 1× 10 11 2× 10 11 5× 10 8 1× 10 9 1× 10 10 3× 10 10 1× 10 11 5× 10 11 1× 10 9 1× 10 10 3× 10 10 1× 10 11 1× 10 12 1× 10 8 1× 10 9 1× 10 10 3× 10 10 1× 10 11 1× 10 12 3×10 8 1×10 9 1×10 10 3×10 10 1×10 11 1×10 12 1×10 13 5×10 3 6×10 3 7×10 3 8×10 3 9×10 3 0 2×10 -8 4×10 -8 6×10 -8 8×10 -8 6x4 8x6 10x8 15x12 40x32 Fig. 9. Trade-off curves between the M2 of LOS displacements and the roughness of the slip distribution, for different p × q-subfault models. Labels give the weight f of roughness in the cost function. R o u g h n e ss Displacement misfit Vol51,2_3,2008 4-03-2009 10:28 Pagina 472 0 1 2 3 4 5 6 7 8 9 10 A lo n g − d ip d is ta n ce (k m ) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 Along−strike distance (km) Best−fit slip distribution, HHS 12x10;f=1x108 0 1 2 3 4 5 6 7 8 9 10A lo n g − d ip d is ta n ce (k m ) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 Along−strike distance (km) 12x10;f=1x1012 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10 11 12 13 Along−strike distance (km) 12x10;f=3x1010 0 1 2 3 4 5 6 7 8 9 10A lo n g − d ip d is ta n ce (k m ) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 6x4;f=3x1010 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10 11 12 13 473 Inversion of synthetic geodetic data for the 1997 Colfiorito events: clues on the effects of layering, assessment of model parameter PDFs, and model selection criteria Fig. 10. Best fit 40 × 32-subfault slip distribution (f = 3 × 1010 gray contours) over the distributed-slip source fault. Fig. 11. Gray contours, all plots: best fit 40 × 32-subfault slip distribution (f = 3 × 1010). Black contours: top left, best fit 12 × 10-subfault slip distribution (f = 3 × 1010); top right, best fit 6 × 4-subfault slip distribution (f = 3 × 1010); bottom left, best fit 12 × 10-subfault slip distribution (f =108); bottom right, best fit 12 × 10-subfault slip distribution (f = 3 × 1012). Vol51,2_3,2008 4-03-2009 10:28 Pagina 473 474 A. Amoruso and L. Crescentini estimation since it is somewhat more conserva- tive in smoothness and computationally less in- tensive than the cross-validation technique. The trade-off curve between the M2 misfit of LOS displacements and the roughness of the slip dis- tribution is shown in fig. 9. It is noticeable that f � 3 × 1010 is in the knee region of the trade- off curve independently of the number of sub- faults (as far as the number of subfaults is large enough to introduce roughness) suggesting that f can be chosen using a single trade-off curve (i.e. a single p x q configuration) for each real case. Figure 10 shows the retrieved 40 × 32 slip distribution overlying the true one. Once again, the most relevant feature consists in the expect- ed focusing effect in the HHS, since the maxi- mum slip patch is shifted by about 1400 m up- ward along the dip direction, i.e. 1000 m in depth. It is important to realize if so many sub- faults (40 × 32) are really necessary to represent the slip distribution or, on the other end, if they are enough. We computed AICc for 6 × 4, 8 × 6, 10 × 8, 12 × 10, 15 × 12, 20 × 16, 40 × 32 subfaults. Minimum AICc is obtained for 12 × 10 subfaults. Because AICc is on a relative scale, here we give the AICc differences with respect to the 12 × 10 case, i.e. 739 for 6 × 4, 84 for 8 × 6, 10 for 10 × 8, 93 for 15 × 12, 384 for 20 × 16, 3778 for 40 × 32. Figure 11, upper plots, shows the capability of AICc to find the model that best explains the data with a mini- mum of parameters, by comparing contours of the slip distribution for 6 × 4, 12 × 10, and 40 × 32, subfaults. Lower plots in fig. 11 show the importance of using the correct weighting fac- tor from trade-off curves when combining dif- ferent terms in the cost function: if f is too small (left plot) the retrieved slip distribution shows unreal small-scale features, if f is too large (right plot) the retrieved slip distribution is too smooth. The same approach should be followed in joint invertion of different kinds of data (see e.g. Amoruso et al., 2008, for EDM and level- ling data). 3. Conclusions Inversion of synthetic coseismic displace- ment data representative of the main event of the 1997 Umbria-Marche sequence has been used to enlight the effects of a crustal layering proper to the Colfiorito area on inferred source features when layering is not taken into ac- count. As a test case we use a recently pub- lished distributed-slip model (Crippa et al., 2006) and we show that the focusing effect of the fault slip distribution can be as large as 1 km. PDFs of model parameters obtained using bootstrapping and NAB are generally in mutual agreement; apparent larger differences relate to the presence of trade-off between selected mod- el parameters (e.g. depth of the fault upper side and along-dip width). PDFs from NAB are usu- ally more scattered than from bootstrapping, especially for subfault slips. This behaviour is probably related to the cost function roughness: on one side NA could be in some trouble in sampling the parameter space efficiently, on the other we use simplex, which could be unable to escape from local minima, for bootstrapping to reduce computing time. This problem deserves attention in the inversion of real data. When inverting for the slip distribution, we use a roughness weight that appears independ- ent of the number of subfaults, thus suggesting that its optimal value can be assessed by a sin- gle trade-off curve (i.e. a single p x q configu- ration) for each real case. The number of sub- faults, that are necessary and sufficient (from the point of view of the information theory) to account for measurements can be obtained from AICc and in this work is as small as 12 × 10. REFERENCES AKAIKE, H. (1974): A new look at the statistical model identification, IEEE Trans. Auto. Control, 19, 716-723. AMORUSO, A., and L. CRESCENTINI (2007): Inversion of lev- elling data: How important is error treatment?, Geo- phys. J. Int., 171, 1352-1362, doi:10.1111/j.1365- 246X.2007.03585.x. AMORUSO, A., L. CRESCENTINI, and R. SCARPA (2002): Source parameters of the 1908 Messina Straits, Italy, earthquake from geodetic and seismic data, J. Geo- phys. Res., 107, 2080, doi:10.1029/2001JB000434. AMORUSO, A., L. CRESCENTINI, and C. FIDANI (2004): Ef- fects of crustal layering on source parameter inversion from coseismic geodetic data, Geophys. J. Int., 159, 353-364, doi:10.1111/j.1365-246X.2004.02389.x. Vol51,2_3,2008 4-03-2009 10:28 Pagina 474 475 Inversion of synthetic geodetic data for the 1997 Colfiorito events: clues on the effects of layering, assessment of model parameter PDFs, and model selection criteria AMORUSO, A., L. CRESCENTINI, and R. SCARPA (2005): Fault- ing geometry for the complex 1980 Campania-Lucania earthquake from levelling data, Geophys. J. Int., 162, 156-168, doi: 10.1111/j.1365-246X.2005.02652.x. AMORUSO, A., L. CRESCENTINI and G. BERRINO (2008): In- version of gravity changes in a layered structure evi- dences magma intrusion during the 1982-1984 Campi Flegrei unrest, Earth Planet Sci. Letters, 272, 181-188, doi: 10.1016/j.epsl.2008.04.040. ÁRNADÓTTIR, T., and P. SEGALL (1994): The 1989 Loma Pri- eta earthquake imaged from inversion of geodetic data, J. Geophys. Res., 99, 21835-21855. BATTAGLIA, M., and P. SEGALL (2004): The interpretation of gravity changes and crustal deformation in active vol- canic areas, Pure appl. geophys., 161, 1453-1467, doi: 10.1007/s00024-004-2514-5. CAPUANO, P., A. ZOLLO, A. EMOLO, S. MARCUCCI, and G. MILANA (2000): Rupture mechanism and source pa- rameters of the Umbria-Marche mainshocks from strong motion data, J. Seismol., 4, 463-478. CINTI, F.R., L. CUCCI, F. MARRA and P. MONTONE (1999): The 1997 Umbria Marche (Italy) earthquake sequence: Relationship between ground deformation and seismo- genic structure, Geophys. Res. Lett., 26, 895-898. CRESCENTINI, L., and A. AMORUSO (2007): Effects of crustal layering on the inversion of deformation and gravity data in volcanic areas: An application to the Campi Flegrei caldera, Italy, Geophys. Res. Lett., 34, L09303, doi:10.1029/2007GL029919. CRIPPA, B., M. CROSETTO, E. BIESCAS, C. TROISE, F. PINGUE, and G. DE NATALE (2006): An advanced slip model for the Umbria-Marche earthquake sequence: coseismic displacements observed by SAR interferometry and model inversion, Geophys. J. Int., 164, 36-45, doi: 10.1111/j.1365-246X.2005.02830.x. CROSETTO, M., C. C. TSCHERNING, B. CRIPPA, and M. CASTILLO (2002): Subsidence monitoring using SAR interferometry: reduction of the atmospheric effects us- ing stochastic filtering, Geophys Res. Lett., 29, 26-29. DE MARTINI, P.M., N.A. PINO, G. VALENSISE, and S. MAZ- ZA (2003): Geodetic and seismologic evidence for slip variability along a blind fault in the Umbria-Marche 1997-1998 earthquakes (central Italy), Geophys. J. Int., 155, 1-11. FUKUSHIMA, Y., V. CAYOL, and P. DURAND (2005): Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise, J. Geophys. Res., 110, doi:10.1029/2004JB003268. HUNSTAD, I., M. ANZIDEI, M. COCCO, P. BALDI, A. GALVANI, and A. PESCI (1999): Modelling coseismic displace- ments during the 1997 Umbria-Marche earthquake (central Italy), Geophys. J. Int., 139, 283-295. HURVICH, C. M., and C.-L. TSAI (1989): Regression and time series model selection in small samples, Biometri- ka, 76, 297-307. INGBER, L. (1993): Simulated annealing: practice versus theory, Math. Comput. Model., 18, 29-57. LAWSON, C.L., and R.J. HANSON (1995): Solving Least Squares Problems, (Soc. Ind. and Appl. Math., Philadelphia, Pa), pp. 337. LUNDGREN, P., and S. STRAMONDO (2002): Slip distribution of the 1997 Umbria-Marche earthquake sequence: joint in- version of GPS and syntetic aperture radar interferometry data, J. Geophys. Res., 107, doi: 10.1029/2000JB000103. MATTHEWS, M.V., and P. SEGALL (1993): Estimation of depth-dependent fault slip from measured surface de- formation with applicationto the 1906 San Francisco earthquake, J. Geophys. Res., 98, 12153-12163. MEGNA, A. (2007): Distribuzione ed evoluzione dello spo- stamento in sistemi di faglie: applicazione a dislo- cazioni di faglia normale in un mezzo bidimensionale eterogeneo nell’Appennino Centro-Settentrionale, PhD Thesis, University of Urbino, Italy. OKADA, Y. (1985): Surface deformation due to shear and tensile faults in a half-space, Bull. seism. Soc. Am., 75, 1135-1154. PINO, N.A., S. MAZZA, and E. BOSCHI (1999): Rupture di- rectivity of the major shocks in the 1997 Umbria- Marche earthquake (central Italy) sequence from re- gional broad band waveforms (1999), Geophys. Res. Lett., 26, 2101-2104. PRESS, W.H., S.A. TEUKOLSKY, W.T. VETTERLING, and B. P. FLANNERY (1992): Numerical Recipes in C: the Art of Scientific Computing, 2nd edn (Cambridge University Press), pp. 994. SALVI, S., et al. (2000): Modeling coseismic displacements resulting from SAR interferometry and GPS measure- ments during the 1997 Umbria-Marche seismic, J. Seism., 4, 479-499. SAMBRIDGE, M. (1999a): Geophysical inversion with a neighbourhood algorithm: I. Searching a parameter space, Geophys. J. Int., 138, 479-494. SAMBRIDGE, M. (1999b): Geophysical inversion with a neighbourhood algorithm: II. Appraising the ensemble, Geophys. J. Int., 138, 727-746. STRAMONDO, S., et al. (1999): The September 26, 1997 Colfiorito, Italy, earthquakes: Modeled coseismic sur- face displacement from SAR interferometry and GPS, Geophys. Res. Lett., 26, 883-886. WANG, R., F. LORENZO MARTIN and F. ROTH (2006): PS- GRN/PSCMP - a new code for calculating co- and post-seismic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation the- ory, Computers & Geosciences, 32, 527-541. ZUCCHINI, W. (2000): An Introduction to Model Selection, Journal of Mathematical Psychology, 44, 41-61, doi:10.1006/jmps.1999.1276. Vol51,2_3,2008 4-03-2009 10:28 Pagina 475