Luciano Scarfì

Horst Langer

Antonio Scaltrito

Seismicity, seismotectonics and crustal velocity structure of the Messina Strait (Italy)

Keywords: Crustal structure, Focal mechanism, Seismicity, Seismotectonics, Tomography

The Messina Strait is the most important structural element interrupting the southernmost part of the Alpine-Apenninic orogenic belt, known as the Calabro-Peloritan Arc. It is being a narrow fan-shaped basin linking the Ionian Sea to the Tyrrhenian Sea. This region is affected by considerable seismic activity which mirrors the geodynamic processes due to the convergence between the African and the Eurasian plates. In the last four centuries, a significant number of disastrous earthquakes originated along the Arc. Among these, the most noteworthy event occurred on December 28, 1908 (known as the Reggio Calabria-Messina earthquake), in the Messina Strait area and caused a large tsunami and more than 100,000 casualties. In this research we focus on the relationships between the general tectonic setting, which characterize the Messina Strait and adjacent areas, seismicity patterns and the crustal structure. We analyzed a data set consisting of more than 300 events occurring in the years from 1999 to 2007, having a magnitude range from 1.0 to 3.8. This data set was exploited in a local earthquake tomography, by carrying out a simultaneous inversion of both the three-dimensional velocity structure and the distribution of seismic foci. We applied the "tomoADD" algorithm, which uses a combination of absolute and differential arrival times and a concept of self-adapting grid geometry, accounting for ray density encountered across the volume. With this method the accuracy of event locations is improved and velocity structure near the source region is resolved in more detail than standard tomography.

Fault plane solutions were obtained for the major and best-recorded earthquakes. The obtained velocity images highlight vertical and lateral heterogeneities that can be associated with structural features striking from NNE-SSW to NE-SW. These results are consistent with important tectonic elements visible at the surface and the pattern delineated by earthquake locations and focal mechanisms.

A c c e p t e d M a n u s c r i p t

Introduction

A hundred years have now passed since the disastrous earthquake on December 28, 1908, one of the largest events ever recorded in the central Mediterranean [START_REF] Boschi | Catalogue of strong Italian earthquakes[END_REF], destroyed Messina, Reggio Calabria and the adjacent areas.

In addition to the damage due to the impact of the seismic waves, much destruction and loss of life was caused by a tsunami which developed as a consequence of the earthquake.

The Messina Strait is without doubt the most important structural element interrupting the southernmost part of the Alpine-Apenninic orogenic belt, known as the Calabro-Peloritan Arc. (Fig. 1). It forms a branch of the Ionian Sea separating SW Calabria from NE Sicily. Its western margin belongs to an escarpment that stretches over 70 km from Messina towards the Mt. Etna area. In W -E direction, the elevation falls from about 1000 a.s.l. to a depth of more than 1000 b.s.l. in the center of the Strait, which means a total elevation difference of 2000 m along a distance of less than 20 km. The considerable orographical differences in NE Sicily, SW Calabria and in the Ionian Sea are the expression of intense tectonic activity. Indeed, from a seismotectonic point of view, the Calabro-Peloritan Arc region is characterized by considerable geodynamic instability [START_REF] Gasparini | Seismotectonics of the Calabrian Arc[END_REF][START_REF] Monaco | Active faulting in the Calabrian Arc and eastern Sicily[END_REF]CPTI work group, 2004). In the last four centuries, a considerable number of violent earthquakes, with epicentral intensities of 10 to 11 degrees MCS (CFTI-Med 4.0), have occurred along the sector of the Arc stretching from the Gulf of S. Eufemia to the Nebrodi Mountains. A significant aspect of this activity is the devastating impact produced on the ground: massive landslides that narrow and sometimes obstruct valleys, the formation of small lakes, the deviation of watercourses, slumps, and the fracturing and liquefaction of soil [START_REF] Bottari | The 1908 Messina Strait earthquake in the regional geostructural frame work[END_REF][START_REF] Murphy | The geomorphological controls on seismically triggered landslides during the 1908 Straits of Messina earthquake, Southern Italy[END_REF][START_REF] Galli | New empirical relationships between magnitude and distance for liquefaction[END_REF][START_REF] Tertulliani | Fenomeni associati al terremoto della Calabria dell'8 settembre 1905[END_REF]. These dramatic phenomena of ground failure are generally considered indicative of the shallow depth of the earthquake source. Indeed, analyses of the macro-seismic fields of the best-documented events favor the hypothesis of a position of the sources -or part of them -close to the surface (see, e.g. [START_REF] Bottari | The 1908 Messina Strait earthquake in the regional geostructural frame work[END_REF]. However, the lack of surface faulting after the largest shocks makes the understanding of the source location as well as the geometry of the fault system a difficult task. We have no ground proof as to whether the important topographical elements around the Messina Strait form true fault systems, perhaps cutting the entire crust and reaching the upper mantle. The identification of such fault systems is further hindered by the great depth of the Ionian Sea (more than 1000 m in the southern part of the Messina Strait) -clearly, there are difficulties in recognising geological formations and tectonic structures underlying the juvenile sediments covering the sea floor at such depths. As a consequence, the dislocation nature of the 1908 Messina Strait earthquake has been a matter of debate (see web site of DISS). [START_REF] Ghisetti | Recent deformations and the seismogenic source in the Messina Strait (Southern Italy)[END_REF] has suggested a graben-like style of faulting for the Messina Strait. For other authors (e.g. [START_REF] Valensise | A 125 kyr-long geological record of seismic source repeatability: the Messina Straits (southern Italy) and the 1908 earthquake (Ms 71/2)[END_REF], the fault causing the 1908 A c c e p t e d M a n u s c r i p t 3 Messina Strait earthquake is seen as a low angle, blindly ending element dipping east, rather than a steeply westward dipping structure parallel to the northeastern coast of Sicily (e.g. [START_REF] Bottari | The 1908 Messina Strait earthquake in the regional geostructural frame work[END_REF].

Research efforts by various institutions have enabled collecting a sizeable set of seismic records, disclosing the possibility of an in-depth investigation of seismotectonic patterns in the area. In particular, the Istituto Nazionale di Geofisica e Vulcanologia has been running a permanent local seismic network from 2001 on; before this time the network was operated by the "Sistema POSEIDON" since the middle of the 1990's. In a recent paper, [START_REF] Barberi | Crustal seismic tomography in the Calabrian Arcregion, south Italy[END_REF] investigated the 3D velocity structure on a fairly wide scale (40 km horizontal grid spacing and 10 km in vertical direction) in an area covering the Southern Tyrrhenian Sea and the Calabrian Arc region. In this investigation, a positive velocity anomaly was found in the Messina Strait at a depth of about 15 km. A further feature of interest here is a negative anomaly stretching from the area north of Mt Etna in ENE direction into the Ionian Sea and the southern tip of the Italian peninsula.

In this research, we focus on the relationship between the general tectonic setting, which characterizes the Messina Strait and adjacent areas, seismicity patterns and the crustal structure. We used a data set consisting of a total of 360 events with a magnitude range from 1.0 to 3.8 occurring in the period from 1999 to 2007. This data set is exploited in a local earthquake tomography (LET), by carrying out a simultaneous inversion of both the three-dimensional velocity structure and the distribution of seismic foci.

In conventional LET (e.g. using the SimulPS12 code; see [START_REF] Evans | User's manual for simulps12 for imaging Vp and Vp/Vs: a derivative of the "Thurber" tomographic inversion simul3 for local earthquakes and explosions[END_REF] the inversion of velocity parameters and earthquake location is carried out on the basis of absolute arrival time readings. These locations are affected by errors collected over the whole ray path, hence the risk that velocity perturbations outside the well illuminated area are mapped into it. Furthermore, the use of absolute arrival times typically requires the definition of station corrections in order to minimize site effects. One way round this dilemma is to use differential travel times, as proposed for instance by [START_REF] Waldhauser | A double-difference earthquake location algorithm: Method and application to North Hayward Fault, California[END_REF]. For earthquake clusters, with foci lying close to each other, travel time errors due to incorrect velocity models in the volume outside the cluster will essentially cancel out. As a consequence of using relative travel times, the use of station corrections becomes obsolete. At the same time, it can be shown (see e.g. [START_REF] Waldhauser | A double-difference earthquake location algorithm: Method and application to North Hayward Fault, California[END_REF][START_REF] Menke | Absolute earthquake locations with differential data[END_REF]) that the use of differential travel times considerably improves the relative locations and provides at least the same quality of absolute hypocenter locations as do standard location techniques. Conventional LET uses a velocity grid defined in Cartesian coordinates and with a regular cell or grid spacing to be defined a priori. The advantage of working with fixed cells lies in its simplicity. For instance, such velocity grids can easily be used in Monte Carlo experiments (see [START_REF] Raffaele | Tomographic inversion of local earthquake data from the Gioia Tauro basin (south-western Calabria, Italy)[END_REF] in order to examine the statistical stability of the results. However, the ray distribution is typically highly uneven due to non-uniform station configuration, the distribution of hypocenters and ray bending. The regular grid geometry makes it difficult to adapt to A c c e p t e d M a n u s c r i p t 4 the unevenness of the rays with the outcome that some cells may have few or no rays, while others may have a very high ray density sampling. Ideally, the inversion grid should be distributed adaptively to match the resolving power of the data (see [START_REF] Zhang | Adaptive mesh seismic tomography based on tetrahedral and Voronoi diagrams: Application to Parkfield, California[END_REF]. Zhang andThurber (2003, 2005) developed computer codes which allow to bypass the problems encountered with absolute hypocenter locations and inflexible velocity grids. Here, we have been using their code tomoADD. Both absolute and relative arrival time readings can be used. The accuracy of the latter can be improved for multiplet events, where relative arrival times can be obtained using correlation and coherency based techniques, such as the crossspectral method. In the adaptive generation of the velocity grid, new grid points are inserted in a tetrahedron of former grid points if a certain threshold value of weighted ray density (DWS) is reached. At the same time, older grid points may be removed if their DWS falls below a limiting value. For more details, see [START_REF] Zhang | Adaptive mesh seismic tomography based on tetrahedral and Voronoi diagrams: Application to Parkfield, California[END_REF].

Finally, our preference for TomoADD was justified from a comparison of the overall RMS travel time residuals obtained with the classical LET method [START_REF] Evans | User's manual for simulps12 for imaging Vp and Vp/Vs: a derivative of the "Thurber" tomographic inversion simul3 for local earthquakes and explosions[END_REF]. TomoADD gave considerably lower residuals (0.13 s) than classical LET using the SimulPS12 code, where an overall residual of 0.29 s was encountered.

Tectonic Setting

The investigated area belongs to the Calabro-Peloritan Arc, which is part of the Apennine-Maghrebian orogenic belt along the Africa-Europe plate boundary. In this sector, convergence processes, from the Neogene-Quaternary on, led to the formation of three distinct tectonic domains: i) the Tyrrhenian back-arc basin (extensional domain), ii) the Alpine-Apenninic orogenic belt (domain of crust shortening) and iii) the Hyblean Foreland (Fig. 1). The Calabro-Peloritan Arc connects the NW-SE-trending southern Apennines with the WSW-striking Maghrebian thrust zones (Fig. 1). Its recent geodynamic evolution has been closely related to the opening of the Tyrrhenian Sea beginning in the middle Miocene, the ESE-ward drift of the Calabro-Peloritan massif and the subduction of Ionian oceanic lithosphere underneath Calabria [START_REF] Finetti | Geophysical study of the Tyrrhenian opening[END_REF]. Several, often contrasting hypotheses, have been proposed for the current geodynamic of the area, such as active subduction beneath the southern Tyrrhenian Sea [START_REF] Barberi | Volcanism of the Southern Tyrrhenian Sea and its geodynamic implications[END_REF], or rifting dominated by slab sinking and passive subduction [START_REF] Finetti | Geophysical study of the Tyrrhenian opening[END_REF][START_REF] Patacca | Tyrrhenian basin and Appenninic arcs: kinematic relations since late Tortonian times[END_REF]. Debate on the most appropriate model is still open. However, the presence of a subducted slab is marked by the occurrence of intermediate and deep earthquakes (between 50 and 400 km) beneath the southern Tyrrhenian Sea, along a NW-dipping Benioff zone [START_REF] Gasparini | Seismotectonics of the Calabrian Arc[END_REF][START_REF] Anderson | The deep seismicity of the Tyrrhenian Sea[END_REF][START_REF] Giardini | The deep seismicity of the Tyrrhenian Sea[END_REF], and by a seismic high-velocity anomaly in the mantle (e.g. Selvaggi and Chiarabba, 1995). The Messina Strait is the most important structural discontinuity cutting the southern part of the Arc, forming a narrow fan-shaped basin linking the Ionian Sea to the Tyrrhenian Sea. Following [START_REF] Ghisetti | Fault parameters in the Messina Strait (southern Italy) and relations with the seismogenic source[END_REF], it is bounded by high angle normal faults with prevailing N-S to NE-SW orientation, and identified as active during Pliocene and Pleistocene times.

Among the major faults, the Larderia-Curcuraci fault system bounds the eastern margin of Sicily. The Reggio Calabria-Calanna-S. Eufemia and the Armo-Delianuova fault systems are situated on the Calabrian side of the Strait. Other NW-SE and E-W striking fault families contribute to a rather complex structural picture, particularly in Calabria (Fig. 1).

Consequently, long-term evolution and the structural setting of the Strait are still a matter of discussion (see among others: [START_REF] Valensise | A 125 kyr-long geological record of seismic source repeatability: the Messina Straits (southern Italy) and the 1908 earthquake (Ms 71/2)[END_REF][START_REF] Tortorici | Recent and active tectonics in the Calabrian Arc (southern Italy)[END_REF].

Seismological, geodetic and structural data demonstrate that the studied area is characterized by active extensional tectonics, at least in the upper part of the crust between depth from surface to 15 km. Indeed, several authors [START_REF] Tortorici | Recent and active tectonics in the Calabrian Arc (southern Italy)[END_REF]D'Agostino e Selvaggi, 2004;[START_REF] Neri | Tectonic stress and seismogenic faulting in the area of the 1908 Messina earthquake, south Italy[END_REF] claim two different styles of deformation accommodating the Europa-Africa convergence in Sicily and in the Calabro-Peloritan Arc. In their opinion the area to the west of the Aeolian-Tindari-Giardini lineament (ATG) is dominated by a compressional regime presumably induced by plates convergence. The area to the east of the lineament is characterized by NW-SE extension that may be related to an Ionian subduction slab rollback.

Data Set and Seismicity

We analyzed seismicity in the area of the Messina Strait by using the data recorded the local network operated by the "Sistema POSEIDON" and subsequently by Istituto Nazionale di Geofisica e Vulcanologia in eastern Sicily and southern Calabria since the 1990's (Fig. 1b). The network consists of three-component stations. Since 2004, stations with analogical telemetry have been upgraded with modern digital VSAT based data transmission and equipped with three-component broadband sensors (40 s). In order to reduce the azimuth gap, the stations of the nearby Aeolian Islands, Mt. Etna and of the National permanent seismic networks are occasionally included as they provide additional onset readings for the most energetic events. All the stations use the same base time, set by GPS time.

In the area extending from 37.81°N to 38.33°N, and 15.35°E to 16.01°E, about 360 earthquakes (Fig. 2) were recorded between 1999 and 2007. In the hypocenter distribution, with coordinates taken from INGV-CT catalogue (see http://www.ct.ingv.it/Sismologia/analisti/default.asp), we recognize a concentration of foci within the uppermost 15-20 km. A similar limitation of focal depth in the Peloritani Mountains has been explained by [START_REF] Langer | Estimation of an optimum Velocity Model in the Peloritani Mountains -Assessment of the variance of model parameters and variability of earthquake locations[END_REF] in terms of a "Brace-Goetze Lithosphere", where rock rheology is assumed to depend on depth and temperature (for more details on this subject see also [START_REF] Stüwe | Geodynamics of the Lithosphere[END_REF]. Greater depth encountered in the south may be related to the presence of the A c c e p t e d M a n u s c r i p t 6 subducted slab, but this must be read with caution as location uncertainties increase due to the unfavorable network configuration. Most of the earthquakes are located within three zones: i) in southwestern Calabria (inland), ii) in the Strait between Messina and Reggio Calabria and iii) in the Ionian Sea, west of Cape S. Alessio. On the other hand, very few events were located in the inland area, between Cape S. Alessio and the northeastern tip of Sicily.

For our investigation, we filtered the original data set according to quality criteria of location. In particular, we selected only well located events, i. e. those with at least 9 observations (P-and S-phases), root-mean-square (RMS) residuals smaller than 0.35 s and horizontal and vertical location errors lower than 2.5 km and 3.5 km, respectively. We also tested the location stability -using the tomoADD code but keeping the velocity parameters fixed -by shifting the trial hypocenters randomly in the space (see [START_REF] Husen | Accurate hypocentre determination in the seismogenic zone of the subducting Nazca plate in northern Chile using a combined on-/offshore network[END_REF]. This helped identify events for which different locations with equivalent travel-time residuals can be found. Note that the presence of events with unstable locations bears the risk of introducing biases as the inversion process may decrease the travel-time residuals by shifting around the hypocenter coordinates instead of adjusting the velocity model parameters properly. In practice, we have compared the locations with unperturbed starting solutions and the locations with starting solutions to which a random perturbation of up to ±8 km was added. We repeated the test five times and considered the realisations where the difference between the solutions were maximum. In doing so, we have a conservative estimate of the stability of the hypocenter locations by removing events with horizontal or vertical location variations greater than 3 km. All tests revealed fairly stable epicenter determinations for almost all the events. In fact, the differences of the results with non-perturbed starting locations and results with randomly perturbed ones was fairly low (1 km or less for 96% of the well locatable events matching the quality criteria mentioned above). Finally, we obtained a data set of 244 well locatable and stable events which were used in the simultaneous inversion of both the velocity structure and the hypocenter location.

3D Imaging

A total of 1785 P and 1163 S absolute arrival times and about 46250 catalogue-derived differential times were jointly inverted using the tomoADD algorithm [START_REF] Zhang | Adaptive mesh seismic tomography based on tetrahedral and Voronoi diagrams: Application to Parkfield, California[END_REF] for hypocentral parameters and velocity structure.

For the given amount of S-wave data available, however, we limit ourselves to discuss the velocity model for P-waves and address to further studies with an enriched data set, when in particular more S-wave arrival times will be available.

We decided to reduce the weight of S-wave arrivals to 20% for fitting the velocity model, whereas their full weight is maintained for the hypocenter location. In doing so, we avoid that S-wave velocities, which are less constrained than those for P-waves, are shifted around artificially reducing the overall residuals. We started the inversion from a regular horizontal grid, with 5x5 km node spacing covering an area of 60x60 km, based on the P and S ray paths of the selected data obtained by the pseudo-bending method (Um and Thurber, 1987) (Fig. 3).

Vertical grid spacing varies between 2 and 6 km and covers a depth range from the surface to 33 km depth, following the 1D reference velocity model [START_REF] Langer | Estimation of an optimum Velocity Model in the Peloritani Mountains -Assessment of the variance of model parameters and variability of earthquake locations[END_REF]. Higher weighting was applied first to the catalogue absolute data in order to obtain a large-scale result. Then, higher weights were assigned to the differential data to refine the event locations and the velocity structure near the source region. The damping values were selected the on the basis of Lcurves of solution and data variances (see [START_REF] Zhang | Adaptive mesh seismic tomography based on tetrahedral and Voronoi diagrams: Application to Parkfield, California[END_REF].

The number and position of inversion mesh nodes change from one iteration to the next as the method automatically adjusts the inversion mesh nodes according to the data distribution. The procedure implemented in TomoADD is based on the construction of tetrahedral and Voronoi diagrams and the use of the Qhull algorithm (see [START_REF] Zhang | Adaptive mesh seismic tomography based on tetrahedral and Voronoi diagrams: Application to Parkfield, California[END_REF]. In this way, the density of the inversion mesh nodes is higher in volumes where more rays are passing. The strategy applied in tomoADD for the design of the mesh is based on the derivative weight sum (DWS), which quantifies the ray density around each model node. DWS is a key parameter for understanding how well the model is resolved [START_REF] Haslinger | 3-D crustal structure from local earthquake tomography around the Gulf of Arta (Ionian region, NWGreece)[END_REF]. Moreover, the adaptive mesh method considers the different ray distribution between the P and S waves. At the end of the iterative process, the inversion mesh nodes for P and S waves are irregularly distributed, accounting for the heterogeneous distribution of the rays (Fig. 4). As a result, ray sampling for each inversion mesh node is more uniform and the tomographic system proves more stable than in a conventional grid with regular spacing (e.g. [START_REF] Zhang | Adaptive mesh seismic tomography based on tetrahedral and Voronoi diagrams: Application to Parkfield, California[END_REF].

In Fig. 5, we compare the distribution of the DWS for a conventional regular mesh created with tomoDD and the one obtained with tomoADD. In the latter, we notice 1379 nodes with non-zero DWS values, all having DWS values greater than 100. The average DWS value is 390 with standard deviation 133. With regular inversion grid (tomoDD), there are 1359 nodes with non-zero DWS value, among which only 36% of them have DWS values greater than 100. Many inversion grid nodes have DWS values less than 10 and the average DWS value is 286. In other words, the adaptive strategy of mesh creation in tomoADD exploits the available information efficiently, whereas the regular grid yields a highly uneven DWS distribution, which becomes a critical issue when only a limited amount of data is available. Again, we notice relatively high P-velocities in the Mts Peloritani and in the Mt. Aspromonte. A negative anomaly is clearly visible in a northeast-southwest striking stripe from the southern part of the Strait to the fault systems of the Calabrian side. At deeper layers (12 and 15 km depth) we notice a slightly modified picture. In fact, in the central and southern parts of the studied area, we notice two zones of low P-wave velocities, trending about ENE-WSW, separated by a ridge with high P-wave velocities. At a depth of 18 km the aforementioned structures are still identifiable, though the picture is blurred by the worsening illumination at this depth. Cross-sections of V P are displayed in Fig. 7.

Stability Analysis with Monte Carlo Tests

Monte-Carlo tests are popular instruments in parameter studies, used to check the stability of forward and inverse modeling. They have the advantage of making very limited a-priori assumptions about possible errors. However, they entail an additional computational burden as the inversions have to be repeated a number of times. This additional effort becomes affordable with increasing computer capacities. During the Monte Carlo tests, we have added uniformly distributed random disturbances to our input data, i.e. the arrival times and the starting velocity model. In particular, uncertainties of the arrival time readings were accounted for by adding a random perturbation of up to ±50 ms for first onset readings and up to ±100 ms for S-phase arrival times. The starting velocity model has been examined in detail by [START_REF] Langer | Estimation of an optimum Velocity Model in the Peloritani Mountains -Assessment of the variance of model parameters and variability of earthquake locations[END_REF] and was found to be the best available one. Nonetheless, as its parameters make up part of the input information we examined the effect of possible uncertainties. We applied random perturbations in the range of ± 0.l km/s both for the P-and S-wave velocities. We repeated the Monte Carlo test 15 times in order to achieve a reasonable significance. Focusing on grid nodes with DWS of at least 100, we found 50 percent of all nodes (i. e. nodes between the 25 % and 75% quartiles) scatter in the range of ± 110 m/s around the results without random perturbation. The standard deviations for the velocity models give a range of ± 130 m/s.

The travel time residuals are the measure of the goodness of fit for the inverted model. A comparison of residuals obtained with and without randomization helps understand whether the assumed uncertainties of input data are realistic or overestimated. In the first case, we expect that travel time residuals obtained with randomization are equivalent to those without, as the accuracy of input data is essentially unaltered. In the second case, we expect that the quality of inverted models worsens significantly. In our case, all inversions carried out with perturbed input data gave slightly higher travel time residuals (about 10% for both absolute and weighted) than the inversion carried out without randomization. We thus conclude that the assumed uncertainties, with respect to arrival time readings and initial guess A c c e p t e d M a n u s c r i p t 9 of the velocity model, are on the safe side as the random perturbation applied to the input data tends to overestimate the true inaccuracy of input parameters. This implies that the expected scatter of the inverted velocity parameters, as far as being caused by the uncertainty of the input parameters, also represents a conservative estimate. On the other hand, we may consider these values a measure of lower threshold for distinguishable anomalies. Velocity differences in the tomographic image which are of the order of 100 to 130 m/s might merely be an effect of an uncertainty of arrival time picking or the parameters of the initial velocity model.

Synthetic Testing

Overall, travel time residuals, DWS and other parameters furnish first-order diagnostics for the assessment of the resolution, whose real significance depends on the individual conditions of the inversion problem and whose understanding needs intuition. We therefore carried out numerical experiments with synthetic models, which give an immediate and straightforward idea of the inversion stability, potentially insignificant or artificial features and the sensitivity with respect to the choice of the starting model. Here, we have used synthetic travel times calculated for typical test models. These data were then inverted using the same starting model, choice of parameters and control values as for the real data.

One type of synthetic test to assess the structural resolving power is the restoring resolution test [START_REF] Zhao | Tomographic imaging of P and S wave velocity structure beneath northeastern Japan[END_REF]. It consists of calculating the synthetic travel times for both the velocity model and hypocenter locations obtained through the inversion. Additionally, noise is added to both travel times and initial guesses. By comparing the results of the inversion of synthetic travel times with those obtained with the real data, we can estimate the restoring capacity of the data set. A major problem of such a restoring resolution test may be local minima in areas of low resolution, leading to the impression of good resolution in actually unresolved areas. The use of a characteristic model avoids this difficulty [START_REF] Haslinger | 3-D crustal structure from local earthquake tomography around the Gulf of Arta (Ionian region, NWGreece)[END_REF]. A characteristic model contains the size and the amplitude of anomalies seen in the inversion results, but with different sign of the amplitude. It is highly unlikely that such a model will represent the same mathematical minimum as the one obtained with the real data. Therefore, areas of low resolution will be indicated by not restoring the characteristic model.

We performed both synthetic tests. Following the strategy of the characteristic model (see [START_REF] Haslinger | 3-D crustal structure from local earthquake tomography around the Gulf of Arta (Ionian region, NWGreece)[END_REF], we designed a velocity model with anomalies at the same position as the ones encountered from the original model but of opposite sign; i. e. positive anomalies become negative, negative become positive. The absolute deviations are maintained. As in our original inversion we started with the optimum 1D model proposed by [START_REF] Langer | Estimation of an optimum Velocity Model in the Peloritani Mountains -Assessment of the variance of model parameters and variability of earthquake locations[END_REF] and applied TomoADD to the synthetic travel times. In the characteristic test we noted that the position and sign of the A c c e p t e d M a n u s c r i p t 10 velocity anomalies are recovered fairly well, at least for nodes where DWS was higher than 100. However, the amount of anomalies is considerably smaller, with a model variance underestimating the true one (Fig. 8a). Travel time residuals for the inverted test model were about 30 ms.

With the restoration test we followed a similar strategy to that applied by [START_REF] Langer | Estimation of an optimum Velocity Model in the Peloritani Mountains -Assessment of the variance of model parameters and variability of earthquake locations[END_REF]. We calculated synthetic travel times for the inverted velocity structure and hypocenter locations. Subsequently, we applied a random perturbation to travel time. Carrying out the inversion without applying random perturbation to the input parameters we achieved a similar travel time residual as with the characteristic model, i. e. 30 ms. Again, the inverted anomalies are similar to the a-priori ones, with respect to position and sign (Fig. 8b). On the other hand, their amount was found smaller than the a-priori values, leading to smaller model variance. With noisy input parameters (perturbations taken from the Monte Carlo experiments, with up to ±50 ms for P-wave onsets, ±100 ms for S-phase picking) the travel time residuals increase to about 50 ms. The residuals reach 100 ms, when an additional systematic error of up to 0.2, representing for instance station corrections, is assumed. The overall characteristic of results is the same as before, in the sense that inverted and a-priori anomalies match both with respect to their position and sign, whereas the inverted model tends to be considerably smoother than the a-priori one.

3D Location

Besides the velocity model, seismic tomography also introduces modifications of the hypocenter location. Compared to the original locations our revised ones show a major degree of clustering (e.g. see the clusters marked with S1 and S2 in Fig. 9). This is no surprise since tomoADD -similar to HypoDD -is based on travel time differences between the events and a major clustering of hypocenters has already been noted in the literature when travel time differences between event clusters are used instead of single absolute arrival time readings (see, e.g, [START_REF] Waldhauser | A double-difference earthquake location algorithm: Method and application to North Hayward Fault, California[END_REF][START_REF] Scarfì | Relocation of microearthquake swarms in the Peloritani mountainsimplications on the interpretation of seismotectonic patterns in NE Sicily, Italy[END_REF]. In our final locations obtained with tomoADD (Fig. 9), we recognize epicenters nicely aligned along the Reggio Calabria-Calanna-S. Eufemia and the Armo-Delianuova fault systems. A further linear alignment of epicenters, striking in a NW direction is found offshore Cape Taormina.

The clusters of events marked as S1 and S2 form two families of multplet events with very similar waveforms, which can be exploited for high precision relative location analysis. We examined the geometry of the two families in more detail using high precision relative location with the master-event technique (e.g. [START_REF] Scarfì | Relocation of microearthquake swarms in the Peloritani mountainsimplications on the interpretation of seismotectonic patterns in NE Sicily, Italy[END_REF] and performed a Principal Component Analysis (PCA) on the covariance matrix of the relative hypocenter locations (see Figs. 10a and10b). For S2, taking the eigenvalues and corresponding eigenvectors, we infer a cluster forming a N 30° E striking planar element with a dip of 70° W. From the square root of the largest eigenvalue and its eigenvector we obtain a A c c e p t e d M a n u s c r i p t 11 hypocenter cloud elongated about 1100 m (68% value, in terms of confidence given by  one standard deviation) in a more or less horizontal direction. Using the same technique, for the S1 cluster, we found a major elongation in a N 20° W striking direction and a dip of about 80° E. We estimate an extension of the data cloud corresponding to about 360 m. In both S1 and S2 clusters, the smallest extension, measured perpendicular the plane spanned by the two major eigenvectors, is of the order of 140 m. This is close to the possible uncertainty of relative location due to intrinsic inaccuracy of travel time determination (see Scarfi et al., 2005).

Fault Plane Solutions and Tectonic Structures.

Fault plane solutions were obtained for the major and best recorded earthquakes occurring in the years 1999-2007. The solutions are based on first arrival P-wave polarities. The nodal lines were identified using the code published by [START_REF] Reasenberg | FPFIT, FPPLOT and FPPAGE: fortran computer progams for calculating and displaying earthquake fault-plane solutions[END_REF]. On the basis of the fault plane solutions and considering patterns inferred from event locations, we may distinguish four seismotectonic regimes: (i) the Messina Strait, (ii) the Mt. Aspromonte, (iii) the Ionian Sea adjacent to the southernmost coast of Calabria, (iv) the offshore region of Cape Taormina and Cape S.

Alessio (see Fig. 11).

The Messina Strait area is characterized by normal faulting events along NE to NNE striking planes, with almost vertical P-axes, and T-axes striking WNW-ESE. Similar patterns were reported by Scarfi et al. (2005) and [START_REF] Giammanco | Possible role of fluid 461 overpressure in the generation of earthquake swarms in active tectonic areas: The case of 462 the Peloritani Mts. (Sicily, Italy)[END_REF] for the nearby Peloritani Mountains area in northeastern Sicily (see also Fig. 11). The orientation of the nodal planes is also consistent with the orientation of the S2 cluster as well as the earlier mentioned epicenter alignments in Calabria. In the Mt. Aspromonte normal faulting as well as strike slip solutions were found, with vertical P-axes or horizontal P-axes striking NNE. In the figure, it can be seen that nodal planes follow trends of tectonic structures visible at the surface. In the Ionian Sea close to the southernmost coast of Calabria the earthquakes show a prevailing normal faulting mechanism with a minor strike slip component, and an extension in NE direction. In the offshore area of Cape Taormina and Cape S. Alessio, there are two events with clear strike slip movement and two normal faults with T-axes oriented WNW-ESE.

Discussion and conclusions

In this paper we have investigated the three dimensional velocity structure in the Messina Strait area. In this context, we re-evaluated the seismicity patterns by relating them to fault plane solutions and tectonic evidences visible at the surface. Seismic tomography was carried out using the recently developed code tomoADD [START_REF] Zhang | Adaptive mesh seismic tomography based on tetrahedral and Voronoi diagrams: Application to Parkfield, California[END_REF], A c c e p t e d M a n u s c r i p t 12 which gave the best arrival time residuals compared to the more classical codes SimulPS12. TomoADD uses absolute arrival times as well arrival time differences among earthquake clusters, and applies a flexible concept of mesh generation. Compared to conventional local earthquake tomography, tomoADD gives at least the same accuracy of earthquake location and guarantees a more homogeneous illumination of the grid cells. The inversion thus makes a more efficient use of the available data set and improves the reliability of results.

We have assessed the stability of our inverted model with respect to uncertainties of the input parameters by applying a Monte Carlo test. We also analyzed the resolution capabilities of the data set, in particular using the DWS distribution and performing synthetic tests. We found that the velocity structure is well resolved at least for nodes where DWS was equal/higher than 100. This condition is met in layers with depth ranging from 6 km to 18 km. From synthetic tests we infer, that the position and sign of velocity anomalies can be recovered fairly well, however, underestimating the variance of the velocity model. In this light our results are understood as a smooth, i. e., conservative estimate of the true heterogeneity of the structure of the crust.

A major point for discussion concerns the structural setting of the Messina Straits and thus the geometry of the seismogenic faults. From the pattern of high-and low-velocity anomalies visible in the layers at different depths as well as in the cross-sections (Figs 6 and 7), we identify a NE-SW striking stripe with negative velocity anomalies. On the Calabrian mainland, this stripe essentially coincides with the zone delimited by the Reggio Calabria -Calanna -S Eufemia (northern border) and Armo -Delianuova faults (southern border, see Fig. 7, cross-sections A-A'', B-B", C-C" and also F-F"). The continuation of this trend in the Ionian Sea can be readily identified in the southernmost profiles, closer to the NE-Sicily coast. Clear positive P-wave velocity anomalies are found in the northern part of the Messina Strait and along the eastern coast of Sicily. The positive anomaly in the Peloritani Mountains is separated from the one of the Strait by a narrow but evident stripe of relatively low velocities. In the deeper layers (12, 15 and 18 km) we notice a ridge of positive anomalies starting in the area of Cape Taormina, then bending in an ENE direction. The ridge crosses the Ionian Sea and is also found in Calabria. For the sake of clarity, figure 12 shows the topography of the 6.5 km/s V P isosurface which summarizes the main features of the obtained velocity structure. In the figure, we also plot the epicenter of the earthquakes located above the isosurface.

Our locations were carried out for a revised data set, focusing on earthquakes satisfying specific quality criteria with respect to the number of available arrival time readings, stability of the starting solution for the hypocenter location, and the goodness of fit (residual). The critical selection of events is necessary in order not to bias the inversion, but also to avoid artifacts, which could be erroneously interpreted as seismotectonic patterns (see [START_REF] Scarfì | High-precision relative locations of two microearthquake clusters in southeastern Sicily, Italy[END_REF][START_REF] Scarfì | Relocation of microearthquake swarms in the Peloritani mountainsimplications on the interpretation of seismotectonic patterns in NE Sicily, Italy[END_REF].

Seismic foci are mainly located on the boundaries of high/low velocity regions (Figs 6,7 and 12), with a characteristic depth between 5 and 15 km. The distribution of well located hypocenters reveals a number of patterns related to A c c e p t e d M a n u s c r i p t 13 structures inferred from other evidences. Among others we note foci accompanying the Armo -Delianuova as well as the Reggio Calabria fault lineaments. Other event locations are found along alignments parallel to the eastern coast of Sicily. High precision relocation of the multiplet family S2 revealed a similar trend, furnishing further evidence for the presence of a NNE striking element. In the southern part of the investigated area, there is an increase of focal depth up to 40 km for the events south of the Calabrian coast (Fig. 9). Here, a rather irregular distribution of hypocenters makes it difficult to identify underlying tectonic structures. On the other hand, the mapping of geological structures is hindered by the sea cover, reaching depths of over 1000 m in the studied area.

Concerning fault plane solutions, earthquakes located in the Messina Strait, the Peloritani Mountains [START_REF] Scarfì | Relocation of microearthquake swarms in the Peloritani mountainsimplications on the interpretation of seismotectonic patterns in NE Sicily, Italy[END_REF] and adjacent areas on the Calabrian land side reveal a normal faulting mechanism along NNE or NE striking planes, which is consistent with the overall regional trends of tectonic motion, i. e. extension in WNW-ESE direction.

Strike slip mechanism was found for two events located offshore Cape Taormina, but again with P-axis striking NNE.

The focal mechanisms of earthquakes falling in the Ionian Sea, south of the Calabrian coast, reveal a change of the stress field as P-axes striking in NW direction prevail.

In conclusion, the various seismological features -3D velocity structure, seismicity patterns, fault plane solutionsconfirm that important tectonic elements visible at the surface reach a considerable depth, at least down to the levels illuminated by earthquake activity. The graben-like structure is well identifiable in the calabrian inland and continues underneath the Messina Strait. The northeastern coastline of Sicily from Taormina to Messina is accompanied by considerable velocity contrasts well identifiable down to a depth of 12 to 15 km. The Messina Strait as well as the adjacent Perloritani Mountains are characterized by extension perpendicular to the NNE striking elements. Both hypocenter locations and focal mechanisms suggest that seismotectonic characteristics in the southern part differ from the picture found in the Messina Strait and adjacent zones. In particular, focal depths of the events increase and the principal direction of major horizontal stress turns from NNE to NW. [START_REF] Ghisetti | Fault parameters in the Messina Strait (southern Italy) and relations with the seismogenic source[END_REF][START_REF] Lentini | Carta geologica della provincia di Messina, Provincia Regionale di Messina[END_REF]. ATG is the abbreviation for Aeolian-Tindari-Giardini faults lineament. Permanent seismic network is also reported: black triangles for Mt. Etna and Aeolian Islands networks and black/white boxes for northeastern Sicily and southern Calabria. Dashed boxes represent the area studied in this paper.

  A c c e p t e d M a n u s c r i p t 7

Fig. 6

 6 Fig.6shows the three dimensional velocity structure for P-wave velocities in the Messina Strait area. We have marked the zones with a reasonable illumination (DWS > 100) with a yellow line. The shallowest layers, at a depth of 6 and 8 km, are characterized by strong positive anomalies in the northern part of the Messina Strait, as well as in the area of the Mts Peloritani and the Mt. Aspromonte. Positive P-wave anomalies are also found more southward along the eastern coast of Sicily, in particular in the offshore area of the Cape S. Alessio. A negative anomaly occurs on the Calabrian side, in an area lying more or less between the Reggio Calabria-Calanna-S. Eufermia and the Armo-Delianuova fault

Fig. 1 .

 1 Figure Captions

Fig. 2 .

 2 Fig. 2. Map view, N-S and E-W cross sections of the studied area with earthquakes located from 1999 to 2007.

Fig. 3 .

 3 Fig.3. 3D sketch of P-wave ray paths traced in the minimum 1D model[START_REF] Langer | Estimation of an optimum Velocity Model in the Peloritani Mountains -Assessment of the variance of model parameters and variability of earthquake locations[END_REF]. Earthquakes and seismic stations are indicated by red circles and blue triangles, respectively.

Fig. 4 .

 4 Fig. 4. Irregular mesh nodes (triangles) for P-(a) and S-waves (b) at the final iteration. At the bottom of the sketches, the projection of the mesh (crosses) and the map of the area (black contour line) are shown.

Fig. 5 .

 5 Fig. 5. DWS value distribution for P waves for the regular inversion grid (a) and irregular inversion mesh (b) at the final iteration (zero values are not shown).

Fig. 6

 6 Fig.6. a) Map of the studied area with the main structural features (from[START_REF] Ghisetti | Fault parameters in the Messina Strait (southern Italy) and relations with the seismogenic source[END_REF][START_REF] Lentini | Carta geologica della provincia di Messina, Provincia Regionale di Messina[END_REF]; b) P-wave velocity model for six representative layers resulting from the 3D inversion. Contour lines are at 0.2 km/s intervals. Red circles represent the relocated earthquakes within half the grid size of the slice. The zones with DWS>100 are circumscribed by yellow contour lines.

Fig. 7 .

 7 Fig. 7. Vertical sections through the P-wave velocity model. The traces of sections are reported in the sketch map (A-A", ... F-F"). Contour lines are at 0.4 km/s intervals. White curves contour the zones with DWS>100. Relocated earthquakes, within ± 4 km from the sections, are plotted as red circles.

Fig. 8 .

 8 Fig. 8. P-wave velocity distribution obtained from the characteristic (a) and restoring (b) resolution tests (see text for further details).

Fig. 9 .

 9 Fig. 9. Final event locations in map and vertical sections. The main fault systems are also shown in the map.

Fig. 10 .

 10 Fig. 10. Relative locations of S1 and S2 clusters by using the master-event-technique. Map view (a), vertical crosssections (b, c) and 3 D sketch (d).

Fig. 11 .

 11 Fig. 11. Focal mechanisms of the major events of the northeastern Sicily and southern Calabria. Their depth are represented together with the fault plane solutions, e. g. Z = 11 (km).

Fig. 12 .

 12 Fig. 12. Topography of the 6.5 km/s P-wave velocity isosurface. Darker colours indicate a greater depth. Contours are at 2.5 km intervals and the numbers indicate the corresponding depths. Relocated earthquakes, above the isosurface, are plotted as red circles.

  Fig_01

Fig_03A

  Fig_03

Fig_08A

  Fig_08

Fig_09A

  Fig_09

Acknowledgments

We wish to thank H. Zhang for allowing us to use the code TomoADD, our colleagues P. De Gori and A. Amato at INGV-CNT for providing data from National seismic network and two anonymous reviewers for carefully reading the manuscript and their constructive criticism. Helpful suggestions by the editor, George Helffrich, are also highly appreciated.


F Riguzzi

G Pietrantonio

R Devoti

S Atzori

M Anzidei

Volcanic unrest of the Colli Albani (central Italy) detected by GPS monitoring test

The Colli Albani volcanic complex, located in central Italy about 15 km SE of Rome, has been dominated by periodic eruptive histories started about 561 ka and ending with the most recent and voluminous activity of the Albano maar (<70 ka) phase. Earthquakes of moderate intensity, gas emissions and significant ground deformations are the recent evidences of a residual activity. We decided to start a monitoring test by installing as first step three GPS permanent stations on the volcanic structure, in sites easily accessible. The analysis of about two years of GPS observations has evidenced a peculiar velocity pattern of the Colli Albani stations with respect to those located nearby, but outside the volcano edifice. With respect to Eurasia, the horizontal velocities are NE directed with magnitudes of 2.2±1.4 mm/yr (RDPI), 3.0±0.8 mm/yr (RMPO) and 3.3±1.2 mm/yr (NEMI). The uplift rates are determined with minor accuracy and range from 3.3 and 6.0 mm/yr. We used a nonlinear inversion algorithm to determine the best-fit parameters for a Mogi spherical source based on the Levenberg-Marquardt least-squares approach. The best fit is obtained with a source at 4.6 km depth beneath the western flank of the volcano and a volume variation of 3.6 •10 -4 km 3 /yr. This result is in agreement with the volume rate retrieved by PS InSAR technique and rather different from the rate inferred from leveling surveys. Consequently, non-linear trends of the hydrothermal system charge cannot be excluded apriori and the continuous GPS monitoring should be considered a priority in assessing the hazard of the Colli Albani.

Introduction

The Colli Albani volcanic complex is located in central Italy, about 15 km SE of Rome, in an area belonging to the potassic and ultrapotassic Roman Magmatic Province, a northwest-trending chain of volcanoes that developed along the Tyrrhenian Sea margin of Italy during middle and late Pleistocene time [START_REF] De Rita | Geological Map of the Alban Hills Volcanic Complex[END_REF]Trigila, 1995). The volcanic history of the Colli Albani is dominated by periodic eruptive histories started about 561 ka and ending with the most recent and voluminous activity of the Albano maar (<70 ka) phase, that cannot be considered completely extinguished [START_REF] Freda | Eruptive history and petrologic evolution of the Albano multiple maar (Alban Hills, Central Italy)[END_REF][START_REF] Funiciello | The Albano maar lake (Colli Albani Volcano, Italy): recent volcanic activity and evidence of pre-Roman Age catastrophic lahar events[END_REF].

Seismic activity and gas emissions are considered the main risk source for the villages located in the Colli Albani area, whereas it is not yet clear which could be the level of hazard due to slow deformations. In fact, this area is characterized by recurrent seismic activity [START_REF] Feuillet | Coupling between earthquake swarms and volcanic unrest at the Alban Hills Volcano (central Italy) modeled through elastic stress transfer[END_REF]Tertulliani and Riguzzi, 1995;[START_REF] Chiarabba | Three-dimensional velocity structure and earthquakes relocation in the Alban Hills volcano, central Italy[END_REF]Amato et al., 1994 and reference therein); temperature and water composition variations [START_REF] Boni | Hydrogeology, fluid chemistry, and thermalism, in The volcano of the Alban Hills[END_REF][START_REF] Calcara | Geochemical monitoring for seismic surveillance, in The volcano of the Alban Hills[END_REF], Carapezza et al., 2008); gas emissions, CO 2 and in minor part H 2 S [START_REF] Carapezza | Geochemistry of the Albano and Nemi crater lakes in the volcanic district of Alban Hills (Rome, Italy)[END_REF][START_REF] Carapezza | Accidental gas emission from shallow pressurized aquifers at Alban Hills volcano (Rome, Italy): Geochemical evidence of magmatic degassing?[END_REF]Tuccimei et al., 2006;[START_REF] Carapezza | Gas hazard assessment in a densely inhabited area of Colli Albani Volcano (Cava dei Selci, Roma)[END_REF][START_REF] Pizzino | Natural gas hazard (CO 2 , 222 Rn) within a quiescent volcanic region and its relations with tectonic; the case of the Ciampino-Marino area, Alban Hills Volcano, Italy[END_REF][START_REF] Chiodini | Carbon dioxide degassing from the Albani Hills volcanic region, Central Italy[END_REF] and reference therein) and by significant ground deformations detected by high precision leveling and PS-InSAR time series. The mean rate of uplift has been estimated from 3 to 7 mm/yr, from PS-InSAR and leveling surveys respectively, on very different time spans [START_REF] Amato | Recent uplift of the Alban Hills volcano (Italy), evidence for magmatic inflation?[END_REF][START_REF] Salvi | Inflation rate of the Colli Albani volcanic complex retrieved by the permanent scatterers SAR interferometry technique[END_REF].

Seismic swarms originate recurrently in the more recent volcanic structures, the freatomagmatics craters of the West side of the Colli Albani, in particular the Albano crater [START_REF] Amato | The 1989-1990 seismic swarm in the Alban Hills volcanic area, central Italy[END_REF]. Seismic tomography analyses recognized the presence of a low velocity volume, identified as magmatic chamber, located about 6 km below this recent area; on the contrary, a high velocity volume was identified under the oldest part of the volcano [START_REF] Chiarabba | Crustal structure, evolution and volcanic unrest of the Alban Hills, Central Italy[END_REF]). The analysis of the teleseismic A c c e p t e d M a n u s c r i p t receiver functions has shown the crustal and uppermost mantle structure beneath the Colli Albani; the Moho is almost flat and located at 23 km, a shallow limestone layer with a variable thickness is present between 4 and 5 km; however, a local complexity of the crustal structure changing from site to site has been also detected [START_REF] Bianchi | Deep structure of the Colli Albani volcanic district (central Italy) from receiver functions analysis[END_REF].

During the last few years, a series of papers on geochronology [START_REF] Freda | Eruptive history and petrologic evolution of the Albano multiple maar (Alban Hills, Central Italy)[END_REF][START_REF] Marra | The Albano maar (Alban Hills volcanic district, Italy): active or dormant volcano?[END_REF], stratigraphyc and physical properties from down hole drilling [START_REF] Mariucci | Looking into a volcanic area: An overview on the 350m scientific drilling at Colli Albani (Rome, Italy)[END_REF], experiments on the magma composition [START_REF] Freda | Magma-carbonate interaction: an experimental study on ultrapotassic rocks from Alban Hills (Central Italy)[END_REF] and eruption products (Sottili et al., 2008) have shed new light on the history and present physical state of the volcanic structure.

Other works (Funiciello et al., 2002;[START_REF] Funiciello | The Albano maar lake (Colli Albani Volcano, Italy): recent volcanic activity and evidence of pre-Roman Age catastrophic lahar events[END_REF][START_REF] Anzidei | The Albano Maar Lake high resolution bathymetry and dissolved CO 2 budget (Colli Albani volcano, Italy): constrains to hazard evaluation[END_REF] highlighted the relevance of high water level variations and catastrophic withdrawal of the Albano maar lake since prehistoric age as possible indicators of sudden variation of CO 2 flow and upwelling of hydrothermal fluids [START_REF] Carapezza | Geochemistry of the Albano and Nemi crater lakes in the volcanic district of Alban Hills (Rome, Italy)[END_REF].

A comprehensive hazard assessment of active and quiescent volcanoes may strongly depend from the morphology of the edifice structure. The landslide hazard evaluation must take into account the presence of sub-aerial and submerged unstable slopes identified in the more recent craters of the Colli Albani [START_REF] Bozzano | Slope dynamics of Lake Albano (Rome, Italy): insights from high resolution bathymetry[END_REF]), whose detection is necessarily based on the knowledge of a detailed digital terrain model (DTM). Therefore, from 2005 to 2007 three different surveys (high precision bathymetric survey of the Albano lake [START_REF] Anzidei | The dark side of the Albano crater lake[END_REF], airborne laser scanning (ALS) and GNSS/RTK survey have been carried out to obtain a high resolution DTM of the area and an estimation of the water volume variations of the Albano lake [START_REF] Pietrantonio | Morphological updating on the basis of integrated DTMs: study on the Albano and Nemi craters[END_REF][START_REF] Riguzzi | Water level and volume estimations of the Albano and Nemi lakes (central Italy)[END_REF].

Most of these researches have been supported by the Department of the Civil Protection through the project DPC115 V3_1, specifically oriented to the definition of potential hazards and crisis levels of the Colli Albani. Considering gas emissions, seismic swarms and ground deformation the most compelling A c c e p t e d M a n u s c r i p t unrest activities of the volcano, one of the priority of the project has been to study the dynamics of the area and their interferences with human activities, including local-scale ground deformation, stress field, slope stability, recent eruptive processes, crater lake evolution, quaternary mass flows.

In this framework, we decided to start a monitoring test at the Colli Albani area by installing as first step three GPS permanent stations in sites easily accessible. With respect to non-permanent networks [START_REF] Betti | GPS sensitivity analysis applied to non-permanent deformation control networks[END_REF][START_REF] Anzidei | Evidence of active crustal deformation of the Colli Albani volcanic area (central Italy) by GPS surveys[END_REF], the continuous monitoring highly improves the accuracy of site velocity estimations by the analysis of long time series of GPS observations [START_REF] Devoti | New GPS constraints on the kinematics of the Apennine subduction[END_REF]. The recent development and densification of GPS permanent networks in Italy, in particular the realization of the RING network (http://ring.gm.ingv.it, [START_REF] Selvaggi | RING working group, La Rete Integrata Nazionale GPS (RING) dell'INGV: una infrastruttura aperta per la ricerca scientifica[END_REF] offers the opportunity to anchor the velocity field of local networks to internationally defined reference systems, currently the ITRF2005 [START_REF] Altamimi | ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters[END_REF], and then to refer local motion in a regional tectonic context.

As underlined in [START_REF] Salvi | Inflation rate of the Colli Albani volcanic complex retrieved by the permanent scatterers SAR interferometry technique[END_REF], unlike other Italian volcanoes, the Colli Albani has not undergone extensive monitoring for long time periods; this gap must be rapidly filled to lay the basis for a more complete assessment of volcanic hazard in such a strongly inhabited and vulnerable area.

Test sites and GPS data processing

We selected the three test sites in locations easily accessible, with electric power facilities: Monteporzio Catone (RMPO), at the Astronomical Observatory; Rocca di Papa (RDPI) at the INGV Observatory, and Nemi (NEMI), at the Padri Verbiti monastery. RMPO is located on the external belt of the Tuscolano-Artemisio caldera, RDPI on the Faete caldera and NEMI on the border separating the ancient Tuscolano-Artemisio caldera and the recent Nemi maar, where together with the Albano and Ariccia maare, the last hydro-magmatic phase took place (Figure 1).

The monument of RMPO is a reinforced concrete column where is immerged a perpendicular steel rod on which the antenna is screwed on. RDPI and NEMI are realized on reinforced concrete buildings A c c e p t e d M a n u s c r i p t with a particular mounting consisting of a short steel pillar (0.65m height) screwed on a steel benchmark permanently and deeply fixed to the top of the building. The antenna is screwed on the top of the pillar and settled in horizontal position through 2 spherical spirit bubbles and 3 screws. Such monument preserves the horizontal and vertical datum when the antenna has to be changed.

RMPO and RDPI are equipped with receivers LEICA GRX1200PRO and LEIAT504 choke-ring antennas (RMPO with radome LEIS, see Figure 1), NEMI is equipped with receiver TRIMBLE 5700 and Zephyr Geodetic antenna TRM41249.00.

We have analyzed the GPS observations at 30 s sampling rates of the Colli Albani test sites in the framework of the processing of all the Italian permanent stations. Consequently, the whole network is divided into 11 clusters of more than 40 stations each one, sharing 11 common anchor sites, i.e. selected sites based on station performance and geographical distribution, subsequently used as core sites for the cluster combination. NEMI, RDPI and RMPO started to work in the middle of 2006, so their data processing covers about 2.5 years.

The data analysis is performed by the Bernese Processing Engine (BPE) of the software Bernese 5.0 [START_REF] Beutler | Bernese GPS Software[END_REF], using the phase double differences as observables. The IGS precise orbits and Earth's orientation parameters are kept fixed and the absolute elevation-dependent phase centre corrections, provided by IGS, are applied. We have estimated each daily solution in a loosely constrained reference frame, close to the rank deficiency condition. Each loosely constrained solution is realized in an intrinsic reference frame, defined by the observations itself, differing from day to day only for rigid network translations, keeping the site inter-distances always well determined. The constraints for the realization of the chosen reference frame are imposed only a posteriori. The daily loosely constrained cluster solutions are then merged into global daily loosely constrained solutions of the whole network applying a classical least squares approach [START_REF] Bianco | Combination of loosely constrained solutions[END_REF].

A c c e p t e d M a n u s c r i p t

We have estimated the velocity field by a purposely designed software (NEVE) that manages the complete stochastic model. The velocity field is estimated by fitting the loosely constrained time series of daily coordinates with the complete covariance matrix, obtaining a loosely constrained velocity solution. We simultaneously estimate site velocities together with annual signals and sporadic offsets at epochs of instrumental changes, according to the following functional model

y A ) sin( 0 ⋅ ≡ ⋅ ∆ + + ⋅ + ⋅ + = i i i i H x t t r x x ϕ ω α
where x 0 is the constant, r the rate, α and φ are respectively the amplitude and phase of the annual signal and H is the Heaviside function useful to detect eventual coordinate offsets (∆x) in the time series.

The unknown parameters of the least squares problem are the components of the vector

T x r x y ) , , , , ( 0 ∆ = ϕ α
and its estimation reads

x C A A) C (A y ˆ1 x T 1 1 x T - - - =
where A is the design matrix, x the observation vector and C x is the covariance matrix of the observations. The loosely constrained velocity field is then minimally constrained and transformed into the ITRF2005 reference frame [START_REF] Altamimi | ITRF2005: A new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters[END_REF] by a 14 parameter Helmert transformation. The errors associated to the velocities derive from the direct propagation of the daily covariance matrices and from the final minimal constraints operator on the loosely constrained velocity field. Concerning the whole network, the most recent sites (2 years of observations) show an uncertainty of 1.5 mm/yr and 3.5 mm/yr respectively in the horizontal and vertical components reaching values down to 0.3 mm/yr and 0.9 mm/yr for the long lasting sites. These errors reflect consistently the expected repeatability of regional velocity field solutions obtained with different methods (e.g. [START_REF] Devoti | New GPS constraints on the kinematics of the Apennine subduction[END_REF][START_REF] Serpelloni | Crustal velocity and strain-rate fields in Italy and surrounding regions: new results from the analysis of permanent and non-permanent GPS network[END_REF]. A more effective representation of the velocity field estimated with this procedure is to show the residual velocities with respect to stable Eurasia (Table 1, Figure 3), whose kinematics is defined through the rotation pole and rate estimated in [START_REF] Devoti | New GPS constraints on the kinematics of the Apennine subduction[END_REF].

The Italian area is the region of the interaction between the African and the Eurasian plates, where the lithosphere is deformed and fragmented in minor blocks, the main of these is the Adriatic plate. The Apennines belt is the chain crossing axially the whole Italian peninsula from NW to SE, it is formed as accretionary prism of the Neogene to present arcuate W-directed subduction of the Adriatic plate. The eastward retreat of the subduction is commonly considered as the most plausible explanation for the progressive eastward migration of the thrust fronts [START_REF] Scrocca | Thrust front segmentation induced by differential slab retreat in the Apennines (Italy)[END_REF]. This retreating is at present well testified by the GPS velocity pattern having the main direction toward NE with respect to the Eurasian reference frame [START_REF] Devoti | New GPS constraints on the kinematics of the Apennine subduction[END_REF]. Then, on average, the Apennine belt constitutes the natural separation of the Italian peninsular velocity field, sites located on the western side display NW-directed velocities whereas sites located on the eastern side display NE-directed velocities.

The Eurasian fixed representation highlights a different trend between the horizontal velocities of the Colli Albani sites (NE-directed) with respect to those located in Rome (NW-directed), about 20 km apart from them (Figure 3).

The three roman sites (INGR, M0SE and ROMA) have a coherent motion with the GPS velocities of sites located along the Tyrrhenian belt, with the above mentioned documented linear trend decreasing from SE to NW [START_REF] Devoti | New GPS constraints on the kinematics of the Apennine subduction[END_REF]. The different behavior of the Colli Albani with respect to the Tyrrhenian sites is better evidenced showing their residual velocities with respect to the roman sites (Table 2, Figure 4).

GPS data modeling

To model the source of deformations, we have referred the GPS velocities with respect to the roman sites (INGR, M0SE, ROMA) with the aim to remove any regional pattern from the local velocity field.

We model GPS data by a Mogi point-pressure source [START_REF] Mogi | Relations between eruptions of various volcanoes and the deformation of the ground surface around them[END_REF]. The underlying reason for the use of such simple model is the low number of available observations (14 measurements from 6 GPS stations, 4 rejected a priori due to clear local effects) that prevents the use of more complex sources like a prolate ellipsoid (Yang et al., 1988) or a sill/dike (Okada, 1985). In fact, though these sources are more realistic and advisable [START_REF] Salvi | Inflation rate of the Colli Albani volcanic complex retrieved by the permanent scatterers SAR interferometry technique[END_REF]Tizzani et al., 2009), they cannot be used with such a low number of observations, unless constraining most of the parameters, either from external observation or from literature.

A Mogi source is defined with only 4 parameters that we can set completely free in the nonlinear inversion. Despite the simplicity of the assumption about the medium (an isotropic and homogeneous half-space with no topography), we can assume that the discrepancy between the Mogi model and a finite element model accounting for the topography falls in the data uncertainty. From [START_REF] Kirchdörfer | Analysis and quasi-static FE modeling of long period impulsive events associated with explosions at Stromboli volcano (Italy)[END_REF] is evident that the last assumption is stronger as the source is deeper: in our modeling the source depth is about 5 km in comparison to highest GPS station located at 700 m.

Under these remarks, we retrieve the source parameters by a non-linear inversion based on the Levenberg-Marquardt (LM) least-squares approach. The LM algorithm is one of the most efficient and widely used optimization algorithm consisting in a combination of a gradient descent and Gauss-Newton iteration [START_REF] Levenberg | A Method for the Solution of Certain Non-Linear Problems in Least Squares[END_REF]Marquardt, 1986). Our implementation of the LM comes from the Fortran MINPACK library (Morè et al., 1980); it is written in IDL and modified with multiple random restarts to guarantee the detection of the global minimum in the optimization process. The cost function is a simple minded weighted mean of the residuals of the form: The best-fit parameters are shown in Table 2, together with the relative standard deviations and the range of values allowed in the inversion. The nonlinear inversion sets the point-source at a depth of about 4.7 km, on the western flank of the Colli Albani complex (see Figure 5). Such results are in a good agreement with those from [START_REF] Salvi | Inflation rate of the Colli Albani volcanic complex retrieved by the permanent scatterers SAR interferometry technique[END_REF], where a prolate ellipsoid double source is used to model DInSAR data.

The displacement is modulated by ∆K, a free parameter that can be converted in terms of volume or pressure variation through the following equations:

π V K ∆ = ∆ µ 3 R P K ⋅ ∆ = ∆
where ∆V is the volume variation, ∆P is the pressure variation for a given value of the radius R. In terms of volume, the best-fit solution corresponds to an increase of 1.14•10 6 m 3 /yr. This solution yields a cost function of 0.62 (the cost function corresponding to the null solution, i.e. the data themselves, is 1.57). The observed and modelled rates are reported in Table 4.

A careful analysis of the uncertainty and the trade-offs affecting the source parameters is also provided. This analysis is performed by perturbing every GPS measurement by a noise source whose value is randomly sampled from a Gaussian distribution with standard deviation equal to that of the considered point. In order to derive how the data uncertainty is mapped to the model parameters, we generate hundreds of perturbed dataset rerunning, each time, the nonlinear inversion. This technique is widely used in the assessment of the model parameters uncertainty [START_REF] Funning | Surface displacements and source parameters of the 2003 Bam, Iran earthquake from Envisat Advanced Synthetic Aperture Radar imagery[END_REF][START_REF] Parsons | The 1994 Sefi Dabeh (eastern Iran) earthquakes revisited: New evidence from satellite radar interferometry and carbonate dating about the growth of an active fold above a blind thrust fault[END_REF]Tizzani et al., 2009) allowing the estimate of the parameter standard deviations and possible trade-offs (Figure 6).

A c c e p t e d M a n u s c r i p t

Discussion and conclusion

We have installed three test sites of continuous GPS monitoring on the Colli Albani volcanic complex to evaluate the velocity components and the achievable accuracies. The sites were selected where power and location facilities made logistically convenient the test. Two observatories and a monastery have satisfied these basic requirements; nevertheless, it has to be underlined this simple network does not assure a suitable azimuthal distribution, since the stations are all located on the eastern flank of the area where the maximum uplift has been detected by PS InSAR technique.

After about 2.5 years of continuous monitoring the results are very encouraging and suggest increasing the number of stations to close the azimuth gap toward SW. The time series analysis of coordinates, after the removal of outliers and instrumental/modeling offsets, has allowed estimating velocities with good accuracy .

With respect to Eurasia, the horizontal velocities of the Colli Albani sites are NE directed with magnitudes of 2.2±1.4 mm/yr (RDPI), 3.0±0.8 mm/yr (RMPO) and 3.3±1.2 mm/yr (NEMI). As known, GPS vertical rates are determined with minor accuracy with respect to the horizontal; in our case we obtain 4.7±5.2 mm/yr (RDPI), 3.3±3.5 mm/yr (RMPO) and 6.0±5.0 mm/yr (NEMI). The accuracies of the rate estimations depend basically on the length of the observation series: long lasting stations can reach accuracies better than mm/yr level, even in the vertical component as shown for INGR and M0SE (Table 1).

The determination of accurate uplift rates on volcanoes is very important, in particular in the Colli Albani area, where vertical motion has been mostly detected till now (Amato and Charabba, 1995;[START_REF] Salvi | Inflation rate of the Colli Albani volcanic complex retrieved by the permanent scatterers SAR interferometry technique[END_REF]. [START_REF] Salvi | Inflation rate of the Colli Albani volcanic complex retrieved by the permanent scatterers SAR interferometry technique[END_REF] have shown that present deformation is concentrated in the area of most recent volcanic activity and the inversion of the observed deformation best fits the data if two From our GPS test, the best-estimated Mogi source (Table 2) is located in the area of the most recent volcanic activity where several phreatomagmatic explosions occurred, generating the Albano, Ariccia and Nemi maare and some coalescent craters. In spite of the poor azimuthal coverage, the estimated parameters (Table 3) are quite in agreement with those of [START_REF] Salvi | Inflation rate of the Colli Albani volcanic complex retrieved by the permanent scatterers SAR interferometry technique[END_REF] and provide indications about the location of future GPS stations.

The horizontal GPS and modeled velocities of the sites located on the volcanic edifice are E-NE directed with magnitudes ranging between 2.7 and 3.3 mm/yr (Table 4, Figure 5). Our data confirm the uplift trend evidenced by other techniques (leveling survey and PS-InSar), but if taken singularly, at this time, appear poorly significant. Nevertheless, plotting the uplift rate of each GPS site vs the distance between the site itself and the center of the estimated Mogi source shown in Figure 5, a clear significant decrease is evidenced. The modeled and GPS uplift trends well agree within the errors, even if the predicted rates appear underestimated with respect to those obtained from GPS (Figure 7).

Another source of uplift could be the flexural-isostatic response to unloading due to the water withdrawal of the Albano lake. Based on the integrated DTM and the recent estimated water level values, we have evaluated about 21.7•10 6 m 3 the water volume loss of the Albano lake from 1993 to 2007, with an average rate of about 1.6•10 6 m 3 /yr [START_REF] Riguzzi | Water level and volume estimations of the Albano and Nemi lakes (central Italy)[END_REF]. We have evaluated that the total unloading effect could produce in this area uplift within 1 mm (Turcotte and Schubert, 1982), less than 0.1mm/yr, and then it can be considered negligible for our purposes.

It has been stressed that both modern seismicity and ground uplifts detected by leveling surveys can be explained by the accumulation of magma of about 94•10 6 m 3 during 43 years, approximately 6 km beneath the zone of seismicity (Feuillet el al., 2004). Hypothesizing a linear trend for the magma A c c e p t e d M a n u s c r i p t accumulation, we can infer a volume change rate of about 21.9•10 -4 km 3 /yr, more than five times the values estimated by GPS and PS-InSAR techniques.

However, non-linear trends of the hydrothermal system charge variations cannot be excluded apriori; our speculation is based on the fact that during different time spans different volume rates of the Mogi sources have been inferred from different techniques (Table 3), so that the continuous GPS monitoring should be considered a priority in assessing the hazard of this volcano.

For the first time, continuous GPS observations have detected deformations of the Colli Albani volcanic structure. The best-fit inversion of GPS velocities is obtained with a Mogi source model, whose parameters are estimated with good accuracy. Our results are in agreement with those retrieved by other techniques, in spite of the poor azimuthal coverage of our sites mainly located on the eastern flank of the area of maximum uplift; at the same time we could expect a velocity W-ward trend of hypothetical GPS sites located on the western flank. Figure 7 shows the circles of equal horizontal velocities at intervals of 0.5 mm/yr. The predicted velocities have a radial pattern from the source reaching the maximum value at about 5 km from the Mogi source where 4.5 mm/yr are expected, thus providing clear indications for future GPS site deployments.

Nevertheless, this work represent a first step and after about three years, the aim of this test seems largely fulfilled; the results achieved till now encourage us to establish new GPS permanent stations primarily on the SW flank of Colli Albani volcano, taking into account the predicted pattern.

Ackowledgments

This work has been partially supported by the Dept. of Civil Protection, Project DPC115_V3 Colli Albani. We wish to thank Alessandra Esposito, Angelo Massucci and Sergio Del Mese for their technical support and all the persons involved in the permanent GPS network maintenance. Some of the figures were made by GMT (Wessel and Smith, 1995). Tuccimei, P., Giordano, G. and Tedeschi, M., CO 2 release variations during the last 2000 years at the Colli Albani volcano (Roma, Italy) from speleothems studies, Earth Planet. Sci. Lett., 243, 449-462, 2006. Turcotte, D.L., and Schubert, G., 1982. Geodynamics, Wiley, New York. Wessel, P., Smith, W. H. F., 1995. The Generic Mapping Tools (GMT) Yang, X., Davis, P.M., and Dieterich, J.H., Deformation from inflation of a dipping finite prolate spheroid in an elastic half-space as a model for volcanic stressing, J. Geophys. Res., 93, 4249-4257, doi: 10.1029/JB093iB05p04249, 1988. A c c e p t e d M a n u s c r i p t 3). A linear decrease of rates is shown with increasing distances.

  A c c e p t e d M a n u s c r i p t

Figure 2

 2 Figure2shows the time series of the N, E, Up coordinates in the ITRF2005 reference frame for the

  ,obs and d i,mod are the observed and modelled displacement of the i-th point and σ i is the standard deviation.

  Mogi sources aligned in a N-S direction at 5 and 7 km depth beneath the western flank of the volcano are considered.

  .,Battaglia, M., Zeni, G., Atzori, S., Berardino, P and R. Lanari, (2009) Uplift and magma intrusion at Long Valley caldera from InSAR and gravity measurements, Geology, 37, 1; 63-66; doi: 10.1130/G25318A.1 Trigila, R. (Editor), The volcano of the Alban Hills, Tipografia SGS,Rome, 1995.

  version 3.0. Technical Reference & Cookbook, SOEST/NOAA.

Figure captions Figure 1 :

 captions1 Figure captions

Figure 2 :

 2 Figure 2: Time series of the North, East and Up components of NEMI, RDPI and RMPO

Figure 3 :

 3 Figure 3: Residual velocity field in central Italy with respect to the Eurasian reference

Figure 4 .

 4 Figure 4. Uncertainty (red bars) and trade-offs (scatter plots) of the source parameters. It

Figure 5 :

 5 Figure 5: Residual GPS velocities (black) with respect to the three roman sites (INGR,

Figure 6 :

 6 Figure 6: Uncertainty (red bars) and trade-offs (scatter plots) of the source parameters. It

Figure 7 :

 7 Figure 7: Observed (red triangles) and predicted (blue open triangles) vertical velocities

Figure 8 :

 8 Figure 8: Circles of predicted equal horizontal velocities at 0.5 mm/yr intervals, their

  Figure(s)

Figure

  Figure(s)

Figure

  Figure(s)

Figure

  Figure(s)

Figure

  Figure(s)

Table 1 :

 1 Residual velocity components, semi-major and semi-minor axes of error ellipses, azimuth of the semi-major axis with respect to Eurasia, at 68% confidence level

	SITE	Lon	Lat	V E	±V E	V N	±V N	Up	±Up	Smax Smin Azim
		° E	° N	mm/yr mm/yr mm/yr mm/yr mm/yr mm/yr	mm	mm	°
	AQUI 13.35 42.37	0.62	0.04	2.34	0.05	0.79	0.12	0.05	0.04	179
	AQUN 13.38 42.34	0.09	1.33	1.91	1.10	3.71	3.11	1.33	1.09	98
	CERT 12.98 41.95	1.63	1.42	2.19	0.88	3.79	1.68	1.42	0.88	92
	CSGP 13.59 42.86	0.94	1.39	2.86	1.58	4.41	4.95	1.59	1.38	-14
	FROS 13.35 41.65	0.22	0.97	1.93	1.29	0.97	3.62	1.29	0.97	177
	GUAR 13.31 41.79	1.73	1.46	0.86	1.25	0.59	3.46	1.47	1.24	102
	INGP 13.32 42.38	1.20	1.12	3.52	1.10	4.94	3.34	1.13	1.09	122
	INGR 12.52 41.83 -1.07	0.60	1.20	0.43	1.39	1.07	0.60	0.43	90
	ITRA 14.00 42.66	0.94	0.54	2.99	0.62	2.52	2.39	0.63	0.54	-3
	LNSS 13.04 42.60 -0.06	0.75	1.65	0.68	-2.85	2.54	0.75	0.68	92
	M0SE 12.49 41.89 -1.06	0.57	1.65	0.45	-0.21	1.33	0.57	0.44	93
	NEMI 12.72 41.72	2.09	1.23	2.60	1.13	6.12	5.04	1.23	1.13	-96
	RDPI	12.71 41.76	1.69	1.49	1.37	1.33	4.73	5.19	1.50	1.31	106
	REFO 12.70 42.96	0.14	0.85	1.76	0.43	4.25	2.88	0.85	0.43	91
	RENO 13.09 42.79	1.00	0.82	2.34	0.95	3.49	3.84	0.95	0.82	180
	REPI	12.00 42.95 -0.03	0.64	1.98	0.59	4.21	2.13	0.64	0.59	94
	RETO 12.41 42.78 -0.56	0.79	1.34	0.59	5.94	3.07	0.79	0.59	90
	RIET	12.86 42.41 -0.29	0.67	1.32	0.79	-0.03	2.44	0.79	0.67	177
	RMPO 12.70 41.81	1.12	0.58	2.79	0.81	3.28	3.46	0.81	0.58	180
	ROMA 12.42 41.91 -0.71	0.86	0.89	1.01	-2.12	3.84	1.01	0.86	174
	RSTO 14.00 42.66	1.24	0.39	2.24	0.34	0.60	0.66	0.39	0.34	91
	TOLF 12.00 42.06 -2.12	0.79	1.84	0.62	1.12	1.47	0.79	0.62	92
	UNOV 12.11 42.72 -0.28	0.49	1.04	0.44	5.47	2.33	0.49	0.44	94
	UNTR 12.67 42.56 -0.03	0.85	1.29	1.52	1.86	1.81	1.52	0.85	180
	VITE	12.12 42.42 -0.43	2.70	1.67	0.66	2.99	3.20	2.70	0.66	90
	VVLO 13.62 41.87 -0.38	0.29	1.79	0.32	5.44	0.81	0.32	0.29	180

Table 2 :

 2 Results of the non-linear inversion for the Mogi source. The standard deviation is calculated from the uncertainty affecting the parameters as shown in Figure4(red bars). The ∆K parameter can expressed in terms either of volume or pressure variation (see the text).

	Parameters	∆K (m 3 /yr)	Depth (km)	East a (km)	North a (km)
	Variation range	[5•10 2 , 3•10 6 ]	[0.00,10.00]	[300.00,310.00]	[4610.00, 4630.00]
	Best fit value	3.64•10 5	4.70	303.35	4620.17
	Standard deviation	6•10 4	0.60	0.53	0.70
	a Coordinates are in UTM-WGS84, Zone 33

Table 3 :

 3 Comparison between modeling parameters of Mogi sources

	Model	observation time span (yr)	Lon (°E)	Lat (°N)	Depth (km)	∆V (km 3 /yr)	RMS
	PS-InSAR M2 North (Salvi et al., 2004)	8	12.665 41.751	4.6	2.0 •10 -4	0.57
	PS-InSAR M2 South (Salvi et al., 2004)	8	12.654 41.666	7.2	4.4 •10 -4	0.57
	levelling S1 (Feuillet et al., 2004)	43	12.687 41.745	4.9	21.9 •10 -4	2.00
	GPS this paper	2	12.635 41.709	4.6	3.6 •10 -4	0.12

Table 4 :

 4 observed and modeled rates, dash indicates data excluded from the nonlinear inversion

				GPS rate components and errors			modeled rates
	Site	V E	V N	Up	±V E	±V N	±Up	V E	V N	Up
		mm/yr mm/yr mm/yr mm/yr mm/yr mm/yr mm/yr mm/yr mm/yr
	INGR	0.0	0.0	0.0	0.8	0.6	1.5	-0.51	0.72	0.24
	M0SE	0.0	0.4	-1.6	0.8	0.6	1.7	-0.22	-	-
	ROMA	0.3	-0.3	-3.5	1.0	1.1	4.0	-0.20	-	-
	NEMI	3.2	1.4	4.7	1.4	1.2	5.0	3.26	0.37	2.17
	RDPI	2.7	0.2	3.3	1.5	1.3	5.2	1.94	1.69	1.39
	RMPO	2.2	1.5	1.9	0.6	0.8	3.5	0.66	1.24	0.51

A c c e p t e d M a n u s c r i p t


Soufiane Belhouideg

Manuel Lagache

Prediction of the effective permeability coefficient in porous media using the finite element method

Keywords: porous media, Finite element method, Permeability, Vermiculite

1

Introduction The aim of this study is to develop a comprehensive Computational Fluid Dynamics (CFD) based model for predicting the permeability coefficient of porous media as exfoliated vermiculite. Vermiculite is a micatype mineral, which is used as a substitute for asbestos for gasket application. A finite element based method is presented for evaluating the effective permeability of such porous solids. The main micro geometry characteristics are determined using scanning electron microscopy observations of vermiculite. Then, the 3-D micro-scale geometries of the porous solids are defined with the ANSYS software. The relation between the microstructure of the porous media and the effective permeability is established. Results of numerical simulations are finally compared with experimental permeability on vermiculite gaskets.

2

Micro geometry definition The estimation of the permeability [1-2] of a porous medium always remains a very difficult task because it is necessary to take into account lot of geometric parameters, as the porosity, the tortuosity, the specific surface… A representation of the 3-D micro-morphologies is essential to evaluate the effect of specific parameters for the optimization of permeability [START_REF] Liu | A simple mathematical model for determining the equivalent permeability of fractured porous media[END_REF]. The direct pore-level numerical simulation of the fluid flow governing equations has became feasible with the advent of high performance of computer and the ability to obtain good geometrical representations of the porous matrix (electronic scan microscopy, tomography). In order to evaluate the main parameters of micro geometry of studied exfoliated vermiculite samples, some electronic microscopy scans were performed. The electronic scan microscopy (SEM) gives some information on the structure of porous media, but it does not generally provide direct threedimensional information. However, the SEM images (Fig. 1) of a sample of vermiculite shows that the layers of vermiculite are shaped like small platelets stacked according an average axis. In order to evaluate the pore network, several SEM observations are performed. The pore network results from the lack of some platelets in this stacking.Informative geometric parameters are needed to describe the complex spatial arrangement of pore. For the present study, spatial distributions and dimensions of pores are partly determined with an Image Processing Software and with direct image analysis of SEM images.

Numerical Model

The proposed method has several distinctive features. First, the 3D micro structures of the porous media are constructed with a parametric program which allow to control the pore sizes, pore orientations, interconnections… Using the FEM software ANSYS and the geometric parameters previously determined, a program generates a pore network with a micro geometry as similar as possible with studied material. Figure 2 illustrated and example of the pore network. The porosities without connection with the pore network are deleted before calculation. Because of the field of sealing application of such vermiculite, the following assumptions are made: the flow is laminar, permanent and isotherm, the volume forces and inertia effects are neglected and the solid matrix is incompressible. The pore network is mesh with 3D Fluid-Thermal tetrahedral element. An optimum mesh size, which produces an accurate solution with a reasonable time, was determined with a preliminary mesh sensitive study. In order to simulate the effective behavior of the macroscopic, adequate boundary conditions are applied on the Representative Volume Element (periodic boundary, pressure gap and no slip boundary condition). Finally, the ANSYS CFD code is used to solve the 3D incompressible continuity and Navier-Stokes equations in the pore-space.

4

Results and discussion The CFD results provide an effective permeability for a porous media with a specific porosity network. Two examples of macroscopic permeability are displayed in Table 1. The first results highlight that the permeability depends strongly of the micro geometry of the porous network. For the same pore volume ratio, the permeability could be very different. Furthermore, this permeability is strongly anisotropic for vermiculite with a pore network build up with pores like a platelet oriented around to an average axis. It can be notice that the effective permeability is of the same order of magnitude with the first experimental results on vermiculite.

Conclusion

In this paper, a FEM based model with the study of SEM images is proposed in order to evaluate the macroscopic permeability of a porous media as the vermiculite. As main benefit, such model could estimate the permeability of a porous media according to the micro geometry. One of the perspectives, is to obtain a more precise description of the porosity network by using systematically efficient Image Processing Software on numerous SEM scan in order to have average information in three directions of the pore network.

Figure 1 -

 1 Figure 1 -image SEM of a sample of the vermiculite

Figure 2 -

 2 Figure 2 -Construction of the 3D model: ( a ) solid part of the REV ( b ) network of the pores ( c ) the connected pores