Giuseppe Pasquale
email: gpasqual@unisannio.it

Raffaella De Matteis

Annalisa Romeo

Rosalba Maresca

Earthquake focal mechanisms and stress inversion in the Irpinia Region (southern Italy)

Keywords: Irpinia Region, Seismicity, Focal mechanisms, Stress inversion

The goal of this study was to estimate the stress field acting in the Irpinia Region, an area of southern Italy that has been struck in the past by destructive earthquakes and that is now characterized by low to moderate seismicity. The dataset are records of 2,352 aftershocks following the last strong event: the 23 November 1980 earthquake (M 6.9). The earthquakes were recorded at seven seismic stations, on average, and have been located using a three-dimensional (3D) P-wave velocity model and a probabilistic, non-linear, global search technique. The use of a 3D velocity model yielded a more stable estimation of take-off angles, a crucial parameter for focal mechanism computation. The earthquake focal mechanisms were computed from the P-wave first-motion polarity data using the FPFIT algorithm. Fault plane solutions show mostly normal component faulting (pure normal fault and normal fault with a strikeslip component). Only some fault plane solutions show strike-slip and reverse faulting. The stress field is estimated using the method proposed by

Introduction

The state of stress within the Earth's crust is of particular interest for geologists and geophysicists as it can provide a better understanding of geodynamic processes [START_REF] Bott | The mechanisms of oblique slip faulting[END_REF][START_REF] Mckenzie | The relation between fault plane solutions for earthquakes and the directions of the principal stresses[END_REF]. The most common methods to determine principal stress orientations are that using breakout data from borehole measurements, which can only be performed at near surface depths [START_REF] Mcgarr | State of stress in the Earth's crust[END_REF], and that using slickenside data. However, in most cases, earthquakes do not show surface ruptures. For this reason, the background seismicity and the aftershocks of a strong earthquake represent the only tools to study the state of stress acting in a region at great depth. Consequently, the use of focal mechanisms to estimate the nature of the stress tensor in the seismogenic zone has been frequently used in the past (i.e., [START_REF] Gephart | An improved method for determining the regional stress tensor using earthquake focal mechanism data: application to the San Fernando earthquake sequence[END_REF][START_REF] Rivera | Stress tensor and fault plane solution for a population of earthquakes[END_REF][START_REF] Michael | Determination of stress from slip data: faults and folds[END_REF]Angelier 1990).

In this study, we have analyzed the seismicity in the Irpinia Region to determine the stress field. It belongs to the most active seismic zones of the Southern Apennines characterized by large destructive earthquakes that occurred both in historical and recent times; the last and strongest of these events was the 23 November 1980 earthquake (M 6.9). This earthquake is the best documented Italian seismic event [START_REF] Giardini | Teleseismic observation of the November 23 1980, Irpinia earthquake[END_REF][START_REF] Pezzo | The 23 November 1980 southern Italy earthquake[END_REF][START_REF] Deschamps | Aftershocks of the Campania-Lucania (Italy) earthquake of the 23 November 1980[END_REF][START_REF] Westaway | The earthquake of 1980 November 23 in Campania-Basilicata (Southern Italy)[END_REF][START_REF] Bernard | The Irpinia (Italy) 1980 earthquake: detailed analysis of a complex normal faulting[END_REF][START_REF] Pantosti | Faulting mechanism and complexity of the November 23, 1980, Campania-Lucania earthquake, inferred from surface observations[END_REF], and was structured into three main rupture episodes (known as those of 0 s, 20 s, and 40 s) that were associated with a surface rupture of about 40 km [START_REF] Amato | Inversion of field data in fault tectonics to obtain the regional stress-III. A new rapid direct inversion method by analytic means[END_REF]. The focal mechanism of this event showed a dominant normal faulting.

The stress field in this region, determined using moderate to strong earthquake focal mechanisms and fault and breakout data, reveals that the Southern Apennines is generally ongoing through a NE-SW extension.

The aim of the present study was to obtain information about the stress field using low magnitude earthquakes. In particular, we used aftershocks of the 23 November 1980 earthquake (M 6.9) to determine the best homogeneous stress tensor, and to compare it with results obtained by other methodologies.

The principal stress orientations can be determined directly from earthquake focal mechanisms through the use of inversion techniques. These techniques yield four parameters: three unit vectors which specify the maximum, minimum, and intermediate compressive principal stress axis orientations and a scalar which describes the relative magnitudes of the principal stresses [START_REF] Gephart | An improved method for determining the regional stress tensor using earthquake focal mechanism data: application to the San Fernando earthquake sequence[END_REF][START_REF] Michael | Determination of stress from slip data: faults and folds[END_REF]Angelier 1990). The accuracy of these inversion techniques depends on the uncertainty of the focal mechanisms and the fault/auxiliary plane ambiguity. In addition to these methods, [START_REF] Rivera | Stress tensor and fault plane solution for a population of earthquakes[END_REF] developed a method whose major advance is to use not previously determined focal mechanisms but the original data that is the polarities of the P arrival and take-off angles for source-station pairs. The focal mechanisms are a by-product of the inverse problem. This method could be useful when the number of the polarities is scarce to determine reliable focal mechanisms.

In this study, we have computed the focal mechanisms of earthquakes located in 3D and 1D velocity models in order to estimate the improvement due to the knowledge of a more realistic model.

We have used the technique developed by Michael (1987a) for inverting focal mechanisms. This technique defines the confidence regions on the quantities obtained through a statistical tool that is known as bootstrap resampling.

Geological and structural setting

Intense hydrocarbon exploration in the Southern Apennines has provided numerous seismic reflection and well data, which have favored the definition of structural models of the Southern Apennines and the reconstruction of its tectonic evolution.

The Apennines are the result of a complex sequence of tectonic events that were associated with the collision between the African and European plate [START_REF] Mckenzie | Active tectonics of the Mediterranean region[END_REF]. They are located between the rather aseismic Tyrrhenian Sea to the West and the compressional Adriatic Sea-Dinarides domain to the East. Over a distance of 100 km or less, there is the change from the active extensional domain of the inner Apennines to areas of active folding and thrusting in the Adriatic Sea [START_REF] Pantosti | Faulting mechanism and complexity of the November 23, 1980, Campania-Lucania earthquake, inferred from surface observations[END_REF].

The compressional tectonics started in the Middle Cretaceous period and lasted until the Oligocene with the continental collision. In the Tortonian age, the rift process started and caused the opening of the Tyrrhenian Basin. Then, following the hinge roll-back of the subducting Adriatic plate, the compressional tectonics migrated towards the East [START_REF] Patacca | Post-Tortonian mountain building in the Apennines: the role of the passive sinking of a relic lithospheric slab. In: Boriani, Bonafede M, Piccado GB, Vai eGB (eds) The lithosphere in Italy[END_REF][START_REF] Patacca | Tyrrhenian basin and Apenninic arc: kinematic relations singe Late Tortonian times[END_REF].

The Southern Apennines is a duplex system [START_REF] Patacca | Tyrrhenian basin and Apenninic arc: kinematic relations singe Late Tortonian times[END_REF]) transported over the flexured southwestern margin of the Apulia foreland. This duplex system consists of a complex architecture of carbonate horsts deriving from the Apulia Carbonate Platform, which is overthrust by rootless nappes. The belt is associated with the Tyrrhenian back-arc basin to the West and with the Bradano foredeep to the East. From Late Tortonian to Lower-Middle Pleistocene, the basin-thrust-belt-foredeep system migrated eastwards, and progressively involved both the basin and carbonate platform paleogeographic domains [START_REF] D'argenio | Schema geologico dell'Appennino Meridionale (Campania e Lucania). Atti del Convegno: Moderne Vedute sulla Geologia dell'Appennino Meridionale[END_REF]) (Fig. 1). Investigations into the Southern Apennines have highlighted the extreme structural complexity of the thrust belt that is due to several factors. The tectonic style is variable since the basin facies terrains were mainly subjected to ductile deformation, whereas the carbonate platform sequences underwent brittle deformation [START_REF] Menardi | Deep structure of the Campania-Lucania arc (Southern Apennines, Italy)[END_REF]. More-over, the kinematic evolution was characterized by out-of-sequence, thrust-propagation processes [START_REF] Roure | Growth processes and melange formation in the Southern Apennines accretionary wedge[END_REF]. In addition, the belt was split up into segments, which underwent different tectonic evolution. In particular, the mountain chain can be subdivided into two major arcs: the NNW-SSE-trending, Molise-Sannio arc to the north, and the WNW-ESE-trending, Campania-Lucania arc to the south, which join in the Irpinia Region.

Finally, starting from the Early-Middle Pleistocene, the axial zone of the chain underwent NE-SW extensional tectonics. This change in the tectonic regime was the cause of extensive volcanism on the Tyrrhenian margin of the chain and the development of large extensional and transtensional structures, which crosscut the pre-existing compressional structures, thus further complicating the internal geometry of the thrust belt. This extensional stress regime causes the seismicity of the Southern Apennines, with large earthquakes centered on the axis of the belt, with depths of mainly around 20 km and characterized by normal faulting focal mechanisms.

Earthquake locations

In this work, we analyzed the earthquakes (M < 4) that occurred immediately after the 23 November 1980 main shock (M 6.9). These were recorded by the National Institute of Geophysics and Volcanology (Istituto Nazionale di Geofisica e Vulcanologia; INGV) seismic network and by a temporary mobile seismic network which was installed in the epicentral area a few days after the mainshock (Fig. 2).

The dataset consists of 2,352 earthquakes recorded by seven stations, on average, resulting in a total of about 21,000 P-wave and 6,700 S-wave arrival time picks and 7,620 P-wave polarity recordings. The location of the earthquakes needs a priori the knowledge of the crustal velocity model. Various one-dimensional (1D) velocity models have been determined in different studies: [START_REF] Deschamps | Aftershocks of the Campania-Lucania (Italy) earthquake of the 23 November 1980[END_REF], [START_REF] Bernard | The Irpinia (Italy) 1980 earthquake: detailed analysis of a complex normal faulting[END_REF][START_REF] Bernard | The Irpinia (Italy) 1980 earthquake: detailed analysis of a complex normal faulting[END_REF][START_REF] Amato | Inversion of field data in fault tectonics to obtain the regional stress-III. A new rapid direct inversion method by analytic means[END_REF] used aftershock arrival times of the 1980 earthquakes to constrain velocities; [START_REF] Improta | Inference on the upper crustal structure of the Southern Apennines (Italy) from seismic refraction investigations and subsurface data[END_REF] determined the velocity structure by interpretation and integration of gravity data, seismic reflection lines, deep wells, and subsurface constraints for the Irpinia Region, up to 7 km in depth. We located the earthquakes using the different 1D velocity models available in literature. The 1D velocity model that minimizes the residuals between the theoretical and observed first P-wave and S-wave travel times is that proposed by [START_REF] Bernard | The Irpinia (Italy) 1980 earthquake: detailed analysis of a complex normal faulting[END_REF].

More recently, [START_REF] Romeo | Crustal structure and stress field for southern Apennines Region from local earthquakes[END_REF] computed a 3D P-wave crustal velocity model in the epicentral area of the 1980 Irpinia earthquake using travel times of aftershocks and background seismicity that occurred from 1988 to 2004 in the Irpinia Region. The volume investigated was 144 × 162 × 30 km 3 , and it was discretized with cell sizes of 9 × 9 × 3 km 3 , with the best resolved zone extending from 0 km to about 20 km in depth. This model overlaps with the volume of the present study.

Considering that: (1) the top of the crustal velocity model corresponds to sea level; and (2) in the area, some stations are located over 1,000 m a.s.l., the present study needed to take into account station corrections. For each station, we identified the outcropping stratigraphic-structural units using geological maps [START_REF] Amanti | Carta Geologica d'Italia Interattiva[END_REF]. Then we associated a correction velocity to each station taking into account the results of [START_REF] Improta | Inference on the upper crustal structure of the Southern Apennines (Italy) from seismic refraction investigations and subsurface data[END_REF][START_REF] Improta | An integrated geophysical investigation of the upper crust in the epicentral area of the 1980, M = 6.9, Irpinia earthquake (Southern Italy)[END_REF]. Finally, on the ground of a vertical ray path, we computed station delays.

The hypocenter parameters were computed again using the [START_REF] Romeo | Crustal structure and stress field for southern Apennines Region from local earthquakes[END_REF] 3D velocity model with the probabilistic, non-linear, global search earthquake location method (Non-LinLoc code; [START_REF] Lomax | Probabilistic earthquake location in 3D and layered models[END_REF]. These authors follow the probabilistic formulation of inversion of [START_REF] Tarantola | Inverse problems = quest for information[END_REF] and [START_REF] Tarantola | Inverse problem theory: methods for data fitting and model parameter estimation[END_REF]. The grid of the probability density function (PDF) values obtained by the grid search represents the complete probabilistic spatial solution of the earthquake location problem. The maximum-likelihood (or minimum-misfit) point of the complete, non-linear location PDF is selected as an "optimal" hypocenter. To make the location program efficient for complex 3D models, the travel times between each station and all of the nodes of an x, y, z spatial grid were computed using a 3D version of the Eikonal finite differences scheme of [START_REF] Podvin | Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools[END_REF]. To compute the travel times, a regular-sized spatial grid of 0.5 × 0.5 × 0.5 km 3 was used.

Figure 2 shows a map of the earthquake locations using the NonLinLoc algorithm and the 3D P-wave velocity model. The average values of location parameters are RMS = 0.4 s, ERH = 3.1 km, and ERZ = 4.3 km. The aftershocks are located along the three fault segments that were associated with the mainshock of the 1980 Irpinia earthquake. Most of the earthquakes occurred within the first 20 km (upper crust).

The earthquake locations using the 3D velocity model do not show significant differences when compared with the locations obtained using the However, the basic difference lies in the computation of take-off angles. Rose diagram in Fig. 3a shows that the take-off angles computed using the 1D layered P-wave velocity model are centered around 30 • and 90 • . This is related to the influence of the layers that make up the velocity model. As for the 1D layered velocity model, also the take-off angles computed using a 1D smooth velocity model (v(z) = (2.27 + 0.4z) km/s) are centered mainly around 30 • and 90 • (Fig. 3b). The take-off angles, computed using the 3D velocity model, are more uniformly distributed (Fig. 3c). We estimated the focal mechanisms for the aftershocks of the 1980 Irpinia earthquake from the P-wave first-motion polarities using the FPFIT algorithm [START_REF] Reasemberg | Fpfit, Fpplot and Fppage. Fortran computer programs for calculating and displaying earthquake fault-plane solutions[END_REF]. FPFIT is a grid search routine that searches for the double-couple fault plane solution that provides the best fit of a given set of first-motion polarities observed for an earthquake.

To compute the high-quality focal mechanisms for the selected dataset of the Irpinia earthquakes, we decided to reject the earthquakes, located in the 3D velocity model, with less than six polarity readings, with horizontal (ERH) and vertical (ERZ) location errors larger than 0.8 km, with an RMS larger than 0.5 s, and with azimuthal gaps larger than 200 • . Worthy of note is that we have considered only the P polarities relative to the P wave arrival time readings with the smaller reading error. This selection reduced the dataset from 2,352 to 321 earthquakes. Moreover, from the computed focal mechanisms, we have selected those mechanisms that had a single solution or at most only one discrepant polarity. For these, among the multiple fault plane solutions relative to the same earthquake, we have selected the focal mechanism which is more consistent with the focal mechanism of the 1980 Irpinia earthquake. In Fig. 4, examples of the quality of the fault plane solutions are shown. The final dataset consists of 139 focal mechanisms (see Table 1 and Fig. 5a).

Mostly, fault plane solutions belong to normal component faulting (pure normal fault and normal fault with a strike-slip component). Only some solutions show strike-slip or reverse faulting. For clearness, the horizontal projections of the extensional (T) axis are shown in Fig. 5b. The length of the small line is inversely proportional to the plunge of the T axis: longer lines represent T axes with smaller values of plunge. A more detailed analysis of these values displayed that For this reason, the earthquake focal mechanisms provide a source of data for stress determination (Michael 1987a). In this work, we used the method devised by Michael (1987a) to determine the stress field from a set of focal mechanisms. Many studies have developed several inverse techniques for finding the stress field orientation using slickenside data (fault orientation and slip direction).

All of these methods are based on two assumptions: first, that the tangential traction on the fault plane should be parallel to the slip direction; and second, that the stress field is uniform within the studied area [START_REF] Gephart | An improved method for determining the regional stress tensor using earthquake focal mechanism data: application to the San Fernando earthquake sequence[END_REF][START_REF] Michael | Determination of stress from slip data: faults and folds[END_REF]Angelier 1990). The use of fault plane solutions is a much more difficult problem with respect to slickenside data because of the fault/auxiliary plane ambiguity. The best algorithm would use a set of fault plane solutions as input and would give the correct choices of fault planes and the best-fit stress tensor as output. [START_REF] Ellsworth | Determination of the stress tensor from focal mechanism data[END_REF] suggested that every possible set of fault plane choices should be considered and inverted. The set of fault planes and slip vectors that fits a single best stress tensor contains the correct choices, and the stress tensor they inferred is the correct stress tensor. This problem was worked out by [START_REF] Gephart | An improved method for determining the regional stress tensor using earthquake focal mechanism data: application to the San Fernando earthquake sequence[END_REF] following a grid search methodology. The algorithm developed by Michael (1987a) is based on that formulation. Applying the Michael technique to a set of focal mechanisms, we obtain the best stress tensor that fits the data and is represented by the orientation of the three principal stress axes, and the scalar which describes the relative magnitudes of the principal stresses and hence constrains the shape of the deviatoric stress ellipsoid, known as parameter φ. This parameter is expressed as:

φ = σ 2 -σ 3 σ 1 -σ 3
where σ 1 , σ 2 , and σ 3 are the maximum, the intermediate, and the minimum compressive principal stress axis, respectively. The basic characteristic of the Michael (1987a) algorithm is the computation of the confidence limits of the principal stress axes directions. Confidence limits are computed by a statistical tool known as bootstrap resampling. From a theoret- we pick data (focal mechanisms) at random from the original dataset. This new dataset will have the same number of data as the original dataset, but will have same mechanisms repeated two or more times while other mechanisms will be absent (Michael 1987a). We inverted this dataset for the stress field, and repeated this process 2,000 times. For example, to define the 95% confidence region, we generated 2,000 new datasets to provide new stress tensors, and the 95% confidence region is Fig. 7 Stereonet plot with the results of the stress inversion using focal mechanism data. The numbers 1, 2, and 3 are the projection of the principal stress axes: σ 1 , σ 2 , and σ 3 , respectively. The thin, dashed, and thick lines circumscribe the 95% confidence regions for σ 1 , σ 2 , and σ 3 , respectively defined by the 95% stress tensors that had the highest normalized scalar product, with the stress tensor determined by the original dataset. The normalized scalar product is a measure of the closeness between two tensors and it is defined as:

3 i=1 3 j=1 B ij C ij 3 i=1 3 j=1 B 2 ij 1 / 2 2 i=1 2 j=1 C 2 ij 1 / 2
where B and C are two tensors.

Empirical tests showed that 2,000 repetitions are enough to yield a stable result (Michael 1987b).

Results

To estimate the stress field in the present study, we analyzed the focal mechanisms of low magnitude earthquakes. Since we cannot decide a priori which nodal plane is the fault plane, we used the Michael (1987a) algorithm, which uses both nodal planes to attempt to choose the correct fault planes while determining the best stress tensor. We also estimated the error on the principal stress axes direction by computing the 95% confidence regions through resampling the data 2,000 times with the bootstrap procedure.

Figure 7 shows a typical representation of the inversion results: the principal stress axes (σ 1 , σ 2 , σ 3 ) with their confidence regions are plotted into the lower hemisphere stereonet. The results obtained by inverting all of the selected focal mechanisms show a nearly horizontal NE-SW minimum compressive stress axis (σ 3 ) and a maximum compressive stress axis (σ 1 ) that is nearly vertical. Although axis σ 3 is horizontal, with a small error on its dip, there is a large uncertainty for its trend, while axis σ 1 shows a large error both for trend and plunge.

To understand if the reason of these large confidence regions is due to the use of a heterogeneous dataset, we used the multi-inverse method proposed by [START_REF] Yamaji | The multiple inverse method: a new technique to separate stresses from heterogeneous fault-slip data[END_REF], which is a numerical technique that is used to separate stresses from heterogeneous fault-slip data. In this way, if our dataset consists of more sets of homogeneous data, the Yamaji method allows us to identify the stress fields compatible with the various subsets. The results are shown in Fig. 8, where the directions of the σ 1 and σ 3 axes are indicated by dots on the lower hemisphere equal-area net to the left and right, respectively. On the left of Fig. 8, the direction of the bars extending from the dots indicates the azimuth of the corresponding σ 3 direction, while its lengths indicate the plunge of the σ 3 axis. In the right panel, the roles of the bars and dots are assigned inversely: the lengths and directions of the bars indicate the σ 1 plunge and azimuth, respectively. Each symbol represents a state of stress and its color represents a stress ratio value (φ). Where the dataset includes two stresses that are compatible with two sets of data, the results show two clusters representing the two stresses. The clusters of dot-bar symbols with the same color and the same bar direction represent significant stresses. In Fig. 8, which was obtained using the selected dataset, the results show a single solution with a large uncertainty, similar to results obtained using the Michael algorithm.

These results show that the dataset is homogeneous and there are no subsets of data that could be the cause of the uncertainty associated with the principal axis directions. We only inverted 115 of the 139 focal mechanisms that show a plunge of the T axis less that 35 • to determine the best stress tensor and we obtain the results shown in Fig. 9. The obtained stress field shows a NE-SW extension, with a horizontal σ 3 (plunge 0 • , trend 230 • ) and a vertical σ 1 (plunge 80 • , trend 320 • ). This stress tensor orientation shows small confidence regions of the principal stresses axes and is consistent with the previous solution obtained using all of the data (see Fig. 7).

These results also show that the large confidence regions of the principal stresses axes, computed using the complete dataset (Fig. 7), are caused only by few data, those with T axis plunge larger than 35 • (24 focal mechanisms). Moreover, the inversion of this small dataset shows an undetermined solution, mainly for the σ 1 and σ 2 axes (Fig. 10). Worthy of note is that the σ 3 axis confidence region is completely not overlapping with that obtained by inverting all the dataset (Fig. 7). This means that the inversion result obtained by inverting only the 24 focal mechanisms is greatly not consistent with that obtained using the complete dataset.

Considering the consistence between the stress field obtained with all of the data and that obtained by inverting the focal mechanisms with a plunge of the T axis ≤35 • , which represent 83% of the complete focal mechanisms, we ascribe the remaining mechanisms to local heterogeneities of the stress field. Figure 11a andb show hypocentral location of earthquakes belonging to both In order to verify if the results obtained using the 3D velocity model are better constrained compared to the results based on 1D velocity models, we have repeated the stress inversion using 1D layered and smooth velocity models. As for the selected earthquakes located in the 3D P waves velocity model, also for events located in the 1D velocity models (layered and smooth), we have computed the focal mechanisms and inverted the selected subsets for stress field inversion. In particular, the results of focal mechanisms inversion of earthquakes located in the 1D smooth velocity model show large uncertainty for orientations of principal stress axes (Fig. 12a) in respect to the results obtained using the earthquakes located in the 3D model. Also, the results obtained using focal mechanisms with plunge of the T axis less than 35 • show large confidence regions, mainly for σ 2 and σ 3 , while σ 1 is better constrained (Fig. 12b). Similar results have been obtained using the 1D layered P waves velocity model (Fig. 12c andd).

The results obtained in this work are compared with those obtained by the use of different methodologies applied recently in Southern Apennines studies and summarized in Fig. 13. In particular, Amato and Montone (1997) and [START_REF] Montone | Active stress map of Italy[END_REF] studied the active stress in [START_REF] Frepoli | Spatial variation in stress in peninsular Italy and Sicily from background seismicity[END_REF]; 5 this study the Southern Apennines through the analysis of breakout and fault plane solution data of moderate to strong earthquakes. They inferred that the Southern Apennines is ongoing through a NE-SW extension, with a horizontal σ 3 and a vertical σ 1 . [START_REF] Frepoli | Spatial variation in stress in peninsular Italy and Sicily from background seismicity[END_REF] used focal mechanisms of the Southern Apennines earthquakes that occurred in the last 20 years to obtain information about the stress field. The inversion of these fault plane solutions, which were mainly normal faults, showed an extensional stress field along a NE-SW direction. Analysis of GPS data also provides a significant SW-NE extension in the Southern Apennines [START_REF] Anzidei | Insights into present-day crustal motion in the central Mediterranean area from GPS surveys[END_REF].

Conclusions

Focal mechanisms of earthquakes from the Irpinia Region have been used to determine the average stress field acting in the area. The dataset is given by aftershocks that follow the 23 November 1980 Irpinia earthquake (M 6.9).

The availability of a 3D velocity model [START_REF] Romeo | Crustal structure and stress field for southern Apennines Region from local earthquakes[END_REF]) allowed us to relocate the earthquakes using a non-linear probabilistic global search algorithm (NonLinLoc-Lomax et al. 2000) to obtain an accurate hypocentral location especially for the finding out the take-off angle values, a basic parameter in the focal mechanisms computation. As a matter of fact, the take-off angles computed using the 3D P-wave velocity model show within a rose diagram a more uniform distribution with respect to those computed using 1D velocity models.

Although the original dataset is large (2,352 earthquakes), our restrictions on the location parameters and the quality of the focal mechanism considerably decreased the amount of data used in this study. Usually, the use of fault plane solutions to determine the stress field shows two different problems: (1) low magnitude earthquakes can be associated with local stress field heterogeneity within a studied area; and (2) the focal mechanism data is characterized by the ambiguity between the fault and the auxiliary planes.

We used the algorithm developed by Michael (1987a), which attempts to choose the correct fault plane while determining the stress tensor, following the assumption that all of the data are consistent with a uniform stress field. The homogeneity of the data was verified using the multiinverse method developed by [START_REF] Yamaji | The multiple inverse method: a new technique to separate stresses from heterogeneous fault-slip data[END_REF]. The large part (115 out of 139) of the focal mechanisms shows a T axis plunge ≤35 • . The results of the inversion show that, using 115 of the 139 focal mechanisms, the best stress tensor obtained is characterized by a horizontal minimum compressive stress axis (σ 3 ) NE-SW oriented and a vertical maximum compressive stress axis (σ 1 ), with small (95%) confidence regions. The remaining data have been explained as focal mechanisms belonging to local heterogeneity of the stress field within the studied area. These data are not consistent with a well constrained stress field, and if they are included in the inversion they increase the uncertainty of the principal stress axes direction.

The comparison of stress tensors inverted on the basis of 1D and 3D velocity models shows that the 3D model represents the main progress providing the better constrained stress field.

The agreement between the results of this study and those obtained using strong to moderate earthquake focal mechanisms and breakout data reveals that low magnitude seismicity represents a useful tool for the determination of stress fields,

Fig. 1

 1 Fig. 1 Simplified geological map of the Southern Apennines showing the major structural features (modified from Tiberti et al. 2005). 1 Molise-Sannio-Lagonegro pelagic basin and related foredeep deposits; 2 Apennine

Fig. 2

 2 Fig. 2 Map of the Irpinia area showing: a aftershock epicenters (circles) and distribution of the stations of the INGV seismic network and the temporary mobile seismic network (triangles). A surface projection of the three fault segments (thick black lines) that ruptured during the 23 November 1980 Irpinia earthquake is also shown (from Pantosti and Valensise 1990). b E-W vertical section showing the hypocenter distribution of the aftershocks

Fig. 3

 3 Fig. 3 Rose diagram for the take-off angles computed using the a 1D layered, b 1D smooth, and c 3D P-wave velocity models. The length of the radius represents the 20% of data

Fig. 4

 4 Fig. 4 Examples of the quality of the fault plane solutions

Fig. 6 a

 6 Fig. 5 a Map showing the distribution of the 139 selected focal mechanisms from the 1980 earthquake epicenter region computed using the P-wave polarity data. The parameters of the computed fault plane solutions are given in Table 1. b Distribution of the focal mechanism T axes. The lengths of the bars are inversely proportional to the plunge of the T axes

Fig. 8 Fig. 11

 811 Fig.8The result of the multi-inverse method[START_REF] Yamaji | The multiple inverse method: a new technique to separate stresses from heterogeneous fault-slip data[END_REF] applied to the data. Each symbol represents a state of stress; the stress ratio is represented by the color, and the direction of the stress axes is indicated by the position and direction of the symbol. The directions of the σ 1 and σ 3 axes are indicated by dots on the lower hemisphere equal-

Fig

  Fig. 12 Stress field results from focal mechanisms of selected earthquakes located in 1D smooth (a) and layered (c) P-wave velocity models. b stress inversion result using the focal mechanisms with T axis plunge ≤35 • for 1D smooth and 1D layered models (d). See Fig. 7 for a description of this figure

Fig. 13

 13 Fig. 13 Stress field orientations obtained from former studies using: 1 GPS data (Anzidei et al. 2001); 2 strong earthquakes (Montone et al. 1999); 3 breakout data (Montone et al. 1999); 4 focal mechanisms of the Southern Apennines earthquakes Frepoli and Amato (2000); 5 this study

Table 1

 1

	Parameters of the 139 selected fault plane solutions	Date 19801201 19801201	T.O. (h/m/s) 02/02/30.45 05/49/4.26	Strike 255 310	Dip 50 30	Rake -60 55	f f
		19801201	23/31/49.03	170	80	165	a
		19801202	04/11/19.68	20	25	-90	f
		19801202	05/16/47.92	215	25	-85	a
		19801202	08/28/32.41	0	40	-110	f
		19801202	21/29/24.36	30	40	-60	f
		19801204	00/41/29.15	30	40	-60	f
		19801204	02/51/31.78	290	50	-5	a
		19801204	06/14/45.35	180	60	-130	a
		19801204	07/12/20.26	290	30	-20	f
		19801204	09/37/17.71	120	80	5	f
		19801204	16/22/5.05	50	75	-100	a
		19801204	18/19/43.47	100	15	85	f
		19801204	20/11/38.28	135	50	55	a
		19801205	02/09/29.51	30	80	-150	a
		19801205	02/38/31.33	280	75	-50	f
		19801205	03/40/12.58	225	20	150	f
		19801205	07/42/31.06	260	80	-105	f
		19801206	02/00/12.3	195	65	-140	a
		19801207	08/44/43.04	325	20	55	f
		19801207	10/08/50.64	255	15	-65	f
		19801207	13/30/5.16	85	70	-65	f
		19801208	04/43/48.78	125	80	-70	f
		19801208	16/06/34.92	335	15	-65	a
		19801208	19/08/46.84	180	55	-90	a
		19801208	20/08/29.26	0	20	-140	f
		19801208	21/21/8.25	320	55	-15	a
		19801208	23/09/58.14	340	60	-160	f
		19801209	02/04/9.51	0	80	-110	a
		19801209	02/31/9.87	95	75	-50	f
		19801209	05/40/33.41	195	35	-175	a
		19801209	09/05/35.2	295	65	20	f
		19801209	14/27/25.52	310	15	-140	a
		19801209	15/22/6.26	85	55	-100	a
		19801209	15/24/29.69	200	50	70	f
		19801209	16/28/4.14	280	30	25	f
		19801209	19/19/4.61	235	15	-115	f
		19801209	23/42/2.82	105	30	5	f
		19801210	01/41/44.28	195	35	20	f
		19801210	02/16/5.94	20	75	-105	a
		19801210	04/50/54.44	230	25	-140	a
		19801210	04/53/14.37	290	45	150	a
		19801210	04/53/58.85	200	70	-75	a
		19801210	07/55/34.51	120	80	20	f
		19801210	12/41/26.56	90	30	-50	f
		19801210	12/57/15	170	40	-110	a
		19801210	14/08/14.38	340	20	180	a
		19801210	14/42/31.3	100	55	-145	a
		19801210	15/34/24.04	190	40	-140	f
		19801210	16/55/32.78	255	20	-120	a
		19801210	17/19/28.64	310	65	25	f
		19801210	19/55/33.79	60	50	-5	a
		19801210	21/18/36.27	15	80	165	a


Francesca Bianco
email: bianco@ov.ingv.it

Lucia Zaccarelli

A reappraisal of shear wave splitting parameters from Italian active volcanic areas through a semiautomatic algorithm

Keywords: Shear wave splitting parameters, Temporal variations, Volcano seismology, Semiautomatic techniques

Shear wave splitting parameters represent a useful tool to detail the stress changes occurring in volcanic environments before impending eruptions. In the present paper, we display the parameter estimates obtained through implementation of a semiautomatic algorithm applied to all useful datasets of the following Italian active volcanic areas: Mt. Vesuvius, Campi Flegrei, and Mt. Etna. Most of these datasets have been the object of several studies (Bianco et al.,

Introduction

Shear wave splitting is the elastic equivalent of the well-known phenomenon of optical birefringence. A shear wave propagating through an anisotropic volume splits into two S-waves (qS 1 and qS 2 ) that travel with different velocities and different polarization directions. This process generates two observables: T d is the time delay between the two split S-waves, and ϕ is the polarization direction of the faster one, qS 1 .

Seismic anisotropy is widely observed in the upper crust regardless of the tectonic setting. This phenomenon has been interpreted as occurring in zones of fluid-filled cracks, microcracks, or preferentially oriented pore spaces. The time evolution of anisotropic distribution of microcracks due to a differential stress, according to the nonlinear anisotropic poroelasticity (APE) model, is explained by fluid migration along pressure gradients between neighboring microcracks and pores [START_REF] Crampin | Modeling the compliance of crustal rock-II. Response to temporal changes before earthquakes[END_REF][START_REF] Zatsepin | Modeling the compliance of crustal rock-I. Response of shear-wave split ting to differential stress[END_REF]. Stress intensity variations cause modifications of the fluid mobility and consequently of the anisotropic characteristics of the medium. In this framework, the shear wave splitting parameters, T d and ϕ, are indicators of the stress field in the upper crust. Time delay variations are related to changes in microcrack density and aspect ratio, while the polarization directions of the two S-waves interchange (90 • flip of ϕ) when the system reaches the overpressurized regime (Crampin et al. 2004 and references therein).

Shear wave splitting parameters have been studied in several volcanic environments such as Hawaii [START_REF] Booth | Shear-wave polarization alignments recorded above the Kaoiki fault zone, Hawaii[END_REF][START_REF] Munson | Crustal shear wave anisotropy in southern Hawaii: spatial and temporal analysis[END_REF] and Long Valley Caldera [START_REF] Savage | Shear wave anisotropy and stress direction in and near Long Valley caldera, California, 1979-1988[END_REF]). Changes in shear wave splitting have been observed during the 1989 and the 2001 eruption at Mt. Etna, Sicily (Bianco et al. 1998b[START_REF] Bianco | Shear wave splitting changes associated with the 2001 volcanic eruption on Mount Etna[END_REF], before the 1996 eruption of Mount Ruapehu, New Zealand [START_REF] Miller | Changes in seismic anisotropy after volcanic eruptions: evidence from Ruapehu[END_REF][START_REF] Gerst | Seismic anisotropy beneath Ruapehu volcano: a possible forecasting tool[END_REF]), before the 1996 October eruption at Vatnajökull, Iceland [START_REF] Volti | A four-year study of shear-wave splitting in Iceland: 2-Temporal changes before earthquakes and volcanic eruptions. In: Nieuwland R (ed) New insights into structural interpretation and modelling[END_REF], and before the M 3.6 earthquake at Mt. Vesuvius, 1999 October (Del [START_REF] Pezzo | Changes in coda decay rate and shear-wave splitting parameters associated with seismic swarms at Mt. Vesuvius, Italy[END_REF].

The main goals of our analysis are: (1) to improve the splitting parameter estimation for monitoring purposes; (2) to define T d and ϕ background values (i.e., out of eruptive periods); and (3) to look for temporal variation of the two observables and its relationship with volcanic activity.

Regarding the first point, we set up a semiautomatic procedure for measuring the splitting parameters in quasi-real time. Regarding points (2) and (3), we performed shear wave splitting measurements on local earthquakes recorded in the following Italian active volcanoes: Vesuvius, Campi Flegrei, and Etna.

Mt. Vesuvius is a stratovolcano located in the Campanian Plain (southern Italy) that experienced various regimes of eruptive activity from effusive to Plinian. Its last eruption occurred in March, 1944 and was a mild effusive event that concluded a period of almost continuous activity dating back to 1631. Since the 1944 eruption ended, the volcano has started a period of quiescence that seems to indicate a closure of its conduit. Up to the time of writing, Vesuvius activity consists of low seismicity and a little degassing at the crater top. We found splitting parameter variations from the background values during two seismic swarms, in 1996 and 1999, in which the largest earthquakes were M 3.4 and M 3.6, respectively.

Campi Flegrei is a resurgent caldera in southern Italy. Its formation has been related to two main explosive eruptions occurred about 38 and 12 ky ago. The most recent activity in 1538 has been the formation of a 200-m-high spatter cone, preceded by a ground uplift of several meters. After this event, the caldera floor started sinking continuously and aseismically. This general trend is occasionally interrupted by faster resurgence episodes accompanied by earthquake swarms. We measured the splitting parameters during the three uplifts: 1982-1984 (180 cm of deformation); 2000 (6 cm of deformation); and 2005-2006 (3 cm of deformation). We found a complex pattern for both splitting parameters in both spatial and temporal domains.

Mt. Etna rises in the NE part of Sicily, Italy. The volcano has experienced many eruptions, and its almost persistent seismicity makes it a natural laboratory for volcano-seismic studies. We calculated the splitting parameters searching for their variations related to the 1988, 1989, and 2001 eruptions. The background value of ϕ is parallel to the direction of the compressive stress field acting on the area, but it shows some variation especially before or during the eruptions. Also the time delays changed, depending on the periods considered and on the eruptive activity.

It has been possible to test the validity of this new semiautomatic method through the comparison of our results with those obtained in previous works for the same datasets. We interpret the observed time variations of the splitting parameters as indicators of complex stress changes, discussing their role as precursors of critical events such as major earthquakes or eruptions in the investigated areas.

Methodology

We first applied a rigorous selection on the seismic records according to the following criteria: (1) clear shear wave onsets with high signal-to-noise ratios: S/N >6 and (2) incidence angles strictly inside the shear wave window (theoretically 35 • ) to ensure no interaction with the free surface. To obtain the signal-to-noise ratio, we measured the noise amplitude before the P-wave onset on all three components. We estimated the two splitting parameters through the implementation of a semiautomatic algorithm. This method minimizes subjectivity and reduces the analysis time required with no loss in accuracy. The semi-automation, in fact, consists of a first step of P and S manual picking, a procedure that needs the precision of an expert eye, and a second step of determination of the two measures: leading polarization and time delay. There exist many ways for estimating the splitting parameters, all of them have different qualities and validity ranges [START_REF] Crampin | A review of techniques for measuring shear-wave splitting above small earthquakes[END_REF]. We have chosen those with few requirements of intervention by the user, even though our analysis system allows easy interaction with screen images to check the results step by step.

We constructed a Matlab algorithm (versatile and user-friendly) that requires as input the previously selected event seismograms, their hypocenter locations, and P and S arrival times. The code evaluates the covariance matrix of a signal portion and diagonalizes it. The horizontal projection of the eigenvector corresponding to the maximum eigenvalue defines the leading polarization direction, ϕ (expressed in degrees from the North). Both N and E components are then rotated to parallel and orthogonal to this azimuth, allowing the distinction between fast and slow waves. Finally, it cross-correlates these two signals whose maximum value occurs at a time lag equal to the qS 1 -qS 2 delay, T d . Note that all results are checked in particle-motion diagrams (hodograms) to avoid possible cycle skipping and other anomalous effects.

Covariance and cross-correlation techniques need a preliminary assignment of start and final points, marking the signal time window of interest. As discussed by [START_REF] Crampin | A review of techniques for measuring shear-wave splitting above small earthquakes[END_REF], the choice of the time window for analysis is crucial for accuracy of the result. In particular, the first splitting parameter has to be evaluated in a signal window containing only the fast S-wave, i.e., where the polarization is roughly linear. Setting this value with no bias is possible only if we already know what the signals look like. Hence, we run the program with a first window guess and looked at the resulting T d distribution. Then, we calibrated the window length on the characteristic time delay found for any particular group of seismic records. In the second parameter estimation, however, the use of cross-correlation function requires signals longer than a period and short enough to not include other phase arrivals later than the S-waves. Our algorithm estimates the S-wave cycle length from the power spectrum calculated over a time window of length [t p ; t p + t s + 0.5 s]. We check the stability of the results trying different window lengths. Summarizing, our algorithm needs the definition of two time windows calibrated on signal characteristics (an example is in Fig. 1). Such semi-automation allows more accurate and objectively reliable results with no loss in terms of data number, since the selection is operated (manually) in the preliminary step of the analysis. Furthermore, it permits us to process more data at the same time producing results in quasi-real time.

We use this procedure to calculate the shear wave splitting parameters for some dataset collected on active volcanic areas in Italy. Most of these data have been already analyzed [START_REF] Bianco | Shear-wave polarization alignment on the eastern flank of Mt. Etna volcano (Sicily, Italy)[END_REF](Bianco et al. , 1998a(Bianco et al. , b, 1999[START_REF] Bianco | The role of hydrothermal fluids in triggering the July-August 2000 seismic swarm at Campi Flegrei (Italy): evidences from seismological and mesostructural data[END_REF][START_REF] Bianco | Shear wave splitting changes associated with the 2001 volcanic eruption on Mount Etna[END_REF][START_REF] Pezzo | Changes in coda decay rate and shear-wave splitting parameters associated with seismic swarms at Mt. Vesuvius, Italy[END_REF]), constituting a useful validation test for the new algorithm, as we applied it exactly on the same number of data previously selected by those authors, following the same selection rules. In order to compare the time delay values with the previous ones and also those of different volcanic areas, we estimated the anisotropy percentage ξ defined following [START_REF] Savage | Shear wave anisotropy and stress direction in and near Long Valley caldera, California, 1979-1988[END_REF]:

ξ = (T d v s /D)×100
, where T d is the time delay, D the hypocenter distance, and v s the S-wave velocity. The averaged ξ represents an index that measures the mean degree of crustal anisotropy over all raypaths from earthquakes to station, and its general meaning is a quantification of both microcrack system density and aspect ratio in the surrounding rocks.

Mt. Vesuvius

Vesuvius is an active volcano, at present in a quiescent stage, located in the southern sector of the Campanian Plain. The last eruption occurred in March 1944 after about 300 years of semi-persistent activity and was an effusive eruption. Its activity began about 300-500 ky ago [START_REF] Santacroce | A general model for the behaviour of the Somma-Vesuvius volcanic complex[END_REF]. The volcanic edifice consists of an older, repeatedly collapsed strato-volcano (Somma), and a more recent cone (Vesuvius). The summit caldera resulted from several collapses related to Plinian (i.e., 79 a.d. "Pompei" eruption; Santacroce 1983) and/or sub-Plinian events (i.e., 1631 eruption; [START_REF] Rosi | Volcanic hazard. In: Santacroce R. (ed) Somma-Vesuvius[END_REF]. The velocity structure beneath Mt. Vesuvius is well known and has been deduced by seismic tomography at different resolution scales [START_REF] Auger | Seismic evidence of an extended magmatic sill under Mt. Vesuvius[END_REF][START_REF] Scarpa | High resolution velocity structure beneath Mt. Vesuvius J Seismol[END_REF]. Moreover, the volcano seismicity pattern is described in many papers, from which we synthesize briefly the main features of earthquake spatiotemporal distribution: the seismicity appears to extend down to 5 km below the central crater, with most of the energy clustered in a volume spanning 2 km in depth, positioned at the border between the limestone basement and the volcanic edifice. The earthquakes are of volcanotectonic type (VT, [START_REF] Chouet | New methods and future trends in seismological volcano monitoring[END_REF], with fault plane orientations showing a highly non-regular spatial pattern. The spectral content of the P-and S-wave trains of the VT events is compatible with stress drops spanning an interval between 1 and 100 bars (100 bars for the highest duration Magnitude, 3.6) and focal dimensions of the order of 100 m, apparently not scaling with the seismic moment (see Scarpa et al. 2002 and references therein). From 1965 to 2000, the seismicity showed an increase both in number and energy of events. The most energetic events occurred in March-April 1996 swarms when the largest number of earthquakes in the last 50 years were recorded (more than 670 events with 1 < M < 3.4; 0 < depth < 6 km b.s.l.) and during the October 1999 seismic swarm that included a duration Magnitude of 3.6, the maximum since at least 1972, but possibly since last eruption in 1944. Overall, the seismicity of Mt. Vesuvius is characterized by a mean rate of 200 events per year.

Data and results

Only the data from two Lennartz 5800 digital 3D stations (BKE and BKN, Fig. 2; 125 s.p.s. with a dynamic range of 120 dB) passed our selection rules and hence have been used for the shear wave splitting analysis. Records from these stations are available from the beginning of 1993 until the end of 2000. Due to the high-dynamic range, the signals are not clipped in the S-wave time window even for the largest magnitude event.

Unfortunately, the stations were not continuously operating during the analyzed period. By applying the procedure described in the previous section and using a v s = 1.8 km/s for the first 2 km, we obtained ξ = 4% for the background seismicity, i.e., for the seismicity excluding the two main crises of 1996 and 1999. In contrast, the same calculation performed on the data collected during the March-April 1996 and October 1999 swarms yield a value of 8% anisotropy. The polarization of the fast S-wave for the background seismicity is NS-oriented at BKN station and approximately EW-oriented at BKE station (Fig. 3). During the two main swarms, however, ϕ strikes NS at BKE, suggesting the occurrence of a 90 • flip mechanism at the station closest to hypocenter locations.

The obtained values of splitting parameters are in agreement with previous works on the same earthquake data (Bianco et al. 1998a[START_REF] Bianco | The seismic crisis at Mt. Vesuvius during 1995 and 1996[END_REF][START_REF] Pezzo | Changes in coda decay rate and shear-wave splitting parameters associated with seismic swarms at Mt. Vesuvius, Italy[END_REF]). The most prominent feature of CF activity is the worldwide-known phenomenon called "bradyseism" consisting of slow subsidence interrupted by phases of noticeably faster ground uplift. Recently, a major bradyseismic crisis occurred in 1982-1984, during which a net uplift of 1.8 m centered on the town of Pozzuoli was accompanied by more than 16,000 earthquakes (M Max = 4.0) mostly located beneath the Pozzuoli-Solfatara area (Fig. 4) at depths between 0 and 4 km b.s.l. [START_REF] Aster | Seismic investigation of Campi Flegrei Caldera. Volcanic Seismology[END_REF]). These earthquakes were recorded by the first digital high-dynamic range mobile seismic network, providing an invaluable data set, still being analyzed, that has allowed a wide researcher community to investigate both the subsurface velocity structure and the dynamics of the seismic source. 1989, 1994, 2000, and 2004. A characteristic of these movements is that seismicity always accompanies the uplift phases, while the subsidence occurs aseismically [START_REF] Saccorotti | The July-August 2000 seismic swarms at Campi Flegrei volcanic complex, Italy[END_REF]). On July-August 2000, the ground deformation reached a maximum value of about 4 cm, and a swarm of low-magnitude VT earthquakes occurred in a period of 3 months. After August 2000, subsidence renewed, and according to geodetic observation, this phase stopped on November 2004 with the onset of a new uplift episode reaching a level of about 4 cm by the end of October 2006 [START_REF] Troise | Renewed ground uplift at Campi Flegrei caldera (Italy): new insight on magmatic processes and forecast[END_REF]). This last phase shows an uplift rate lower than that of 2000, and at the present days, the very last ground deformation measurements seem to indicate a general cessation of the phenomenon. Remarkable swarms of both VT and LP (Long Period, [START_REF] Chouet | New methods and future trends in seismological volcano monitoring[END_REF] earthquakes have accompanied this last ground uplift, with the highest number of seismic events ever recorded since 1985 [START_REF] Saccorotti | Seismicity associated with the 2004-2006 renewed ground uplift at Campi Flegrei caldera, Italy[END_REF]).

Data and results

We measured the splitting parameters for earthquakes which occurred during the 1982-1984 (dataset 1), 2000 (dataset 2), and 2004-2006 (dataset 3) crises.

Concerning dataset 1, 80 earthquakes recorded at three stations (W11, W12, and W21, see Fig. 4) passed our selection rules; the depths of the selected seismicity mainly spans the first 3-4 km of the upper crust. By applying the semiautomatic algorithm and using a v s = 1.7 km/s for the first 2 km, we calculated that ξ is 7%. Spatially, the polarization of the fast wave shows an irregular pattern: the polarization is approximately NNW/SSE-oriented at W11 and W12 stations, and NNE/SSW-oriented at W21 station (Fig. 4). [START_REF] Savage | Shearwave anisotropy of active tectonic regions via automated S-wave polarization analysis[END_REF] studied the shear wave splitting at Campi Flegrei for 51 microearthquakes that occurred during the 1984 bradyseismic crisis. They too processed data from station W11, located just at the center of the caldera, finding a minimum anisotropy of 7% and a fast shear wave orientation striking NNE-SSW, values that are compatible with our results.

For dataset 2, 19 events recorded at the digital station SLF (Fig. 4) passed the selection rules. By applying the semiautomatic algorithm, we determined that the fast polarization strikes mainly NNW-SSE, confirming the orientation retrieved by [START_REF] Bianco | The role of hydrothermal fluids in triggering the July-August 2000 seismic swarm at Campi Flegrei (Italy): evidences from seismological and mesostructural data[END_REF]. With the same S-velocity of 1.7 km/s for the first 2 km, we calculated ξ to be about 4%. This value is perfectly comparable with those observed at the Californian Geysers geothermal field [START_REF] Evans | Shearwave splitting from local earthquakes at The Geysers geothermal field, California[END_REF] and at Long Valley caldera [START_REF] Savage | Shear wave anisotropy and stress direction in and near Long Valley caldera, California, 1979-1988[END_REF]), but is lower than the value found by [START_REF] Savage | Shearwave anisotropy of active tectonic regions via automated S-wave polarization analysis[END_REF] at Campi Flegrei and in this study for dataset 1.

For dataset 3, 20 earthquakes recorded at three stations (ASB2, AMS2, and BGNB, Fig. 4) passed our selection rules. The three stations are digital high-dynamic Lennartz M-24/Lennartz MarsLite equipped with broadband sensors (LE-3D 20 s/ Guralp CMG 40 s). By applying the semiautomatic algorithm, we determined that the fast polarization strikes approximately NNW-SSE at ASB2 and NE-SW at AMSE and BGNB stations. Using a v s = 1.7 km/s for the first 2 km, we found an anisotropy percentage of 3%, which is the lowest value we retrieved in the Campi Flegrei area.

The difficulty in discussing these results is that we cannot distinguish between background and crisis seismicity, since at Campi Flegrei, as highlighted above, earthquakes occur only during the uplift crisis.

Mt. Etna

Mt. Etna volcanic activity, mainly characterized by effusive basaltic eruptions, began around 0.7 My ago, continuing almost uninterruptedly until the present. The volcano has developed at the intersection of two regional tectonic lineaments, the NNW-SSE trending Hybleo-Maltese escarpment, which separates the thick inland continental crust of the African platform from the Ionian Mesozoic oceanic crust, and the NE-SW Messina-Fiumefreddo fault that marks a rift zone between south Calabria and northeastern Sicily extending as far as the volcanic area. All these tectonic features affect the morphology, especially of the eastern side of the volcano (Bonforte and Puglisi 2006 and references therein). The volcanic activity is always accompanied by the occurrence of seismicity. Earthquakes may also happen independently from magma movements due to the dynamics of the several active faults through the volcano.

Data and results

Shear wave splitting parameters have been investigated several times during different periods of activity. [START_REF] Bianco | Shear-wave polarization alignment on the eastern flank of Mt. Etna volcano (Sicily, Italy)[END_REF] measured the splitting parameters on a data set collected in spring-summer 1988, a period with no eruption. Eleven earthquakes recorded at 10 stations, all located on the eastern sector, passed the selection rules, and the selected events occurred in the first 10-15 km of the upper crust, with magnitude between 2.3 and 3.1. The stations were digital high-dynamic Lennartz PCM 5800 equipped with 1 Hz Mark L4C-3D and 2 Hz Mark L4A-3D sensors. The authors found that the fast polarization at all stations strikes, approximately, in the NS direction, which is the direction of the maximum compressive stress deduced for the area. We applied our semiautomatic procedure on subset of the data investigated by [START_REF] Bianco | Shear-wave polarization alignment on the eastern flank of Mt. Etna volcano (Sicily, Italy)[END_REF], consisting of 10 events recorded at two stations (DMT and POM, Fig. 5). The fast polarization strikes in NS direction at both stations (Fig. 6), as found previously. Using a v s = 2.2 km/s for the first 8 km, we calculated that the background ξ averages 1.5%.

A similar seismic network was installed again on the eastern sector of the volcano during May-November 1989; this period includes the occurrence of the September-October eruption. We calculated the splitting parameters for the signals collected at two stations of this network (NOC and CCV, in Fig. 5). Ten earthquakes passed the selection rules; again, the selected events sampled the first 10-15 km of the upper crust, with magnitude between 2.1 and 3.5. At all the stations, ϕ is confirmed in the NS direction (Fig. 7), while the anisotropy, considering the same mean velocity, averages 3%.

During the 1989 eruption, a NNW/SSEoriented fracture system opened at 2,700 m above sea level (September 27) and propagated 6.5 km downslope to 1510 m a.s.l. (October 3). Bianco et al. (1998b) performed shear wave splitting analyses during the propagation of the fracture system in order to study the correlation between seismic anisotropy and the syn-eruptive cracking process. Due to the selection rules, the available dataset was reduced to 13 earthquakes collected at one (B92) digital high-dynamic Lennartz PCM 5800 equipped with 1 Hz Mark L4C sensor. The qS 1 polarization strikes NE-SW, i.e., almost orthogonal to the strike of the NNW-SSE fractures nearest B92 and subparallel to the ENE-WSW cracks associated to the folds. Interestingly, the delay between the split waves exhibited a time-related trend with the fracture propagation (Bianco et al. 1998b). Applying our semiautomatic procedure, we confirm the ENE-WSW orientation for the fast S-wave (Fig. 7). Using a v s = 1.2 km/s for the first 2 km, the anisotropy for this dataset is around 6%.

Lastly, we applied the semiautomatic algorithm to infer the splitting parameters for 134 earthquakes which occurred during 2001 in a period including the July eruption. This dataset was collected at two digital high-gain 3C seismic stations MNT and ESP (Fig. 5). Our procedure returns for the complete set of data qS 1 polarization directions striking approximately NNW-SSE at both ESP and MNT stations (Fig. 8), in agreement with the results of [START_REF] Bianco | Shear wave splitting changes associated with the 2001 volcanic eruption on Mount Etna[END_REF]. The nearly EW polarizations at MNT occur mostly in the period immediately before the eruption when there are 90 • flips in shear wave polarizations (see Fig. 3b in [START_REF] Bianco | Shear wave splitting changes associated with the 2001 volcanic eruption on Mount Etna[END_REF]). We took a v s = 2.2 km/s over the first 5 km to calculate a ξ value of 4.4%.

The temporal variation of the splitting parameters during the occurrence of the 2001 eruption has been widely investigated, for the same dataset, by [START_REF] Bianco | Shear wave splitting changes associated with the 2001 volcanic eruption on Mount Etna[END_REF]. These authors found that the normalized time delays show a systematic increase starting several days before the onset of the eruption (July 17). The time delays also exhibit a sudden decrease shortly before the start of the eruption, suggesting the possible occurrence of a 'relaxation' phenomenon due to crack coalescence. This behavior has many similarities to that observed before a number of earthquakes elsewhere [START_REF] Gao | Observations of stress relaxation before earthquakes[END_REF]). On several occasions before Etna activity onset, the qS 1 polarization directions at the stations closest to the eruption show 90 • flips, meaning that faster and slower shear waves exchange their polarizations. In particular, the unflipped ϕ strikes approximately in NS direction, as always observed at Mt. Etna.

Discussions and conclusions

Three Italian active volcanoes have been investigated in order to assess their anisotropic features through the implementation of a semiautomatic algorithm applied to all the available datasets and also to check the association between temporal variations of splitting parameters and the occurrence of changes in the dynamic pattern of the volcanic activity. As recommended by [START_REF] Gao | SWAS: a shear-wave analysis system for semi-automatic measurement of shear-wave splitting above small earthquakes[END_REF], the semiautomatic techniques need user control at a particular step of the analysis (in our case the window length calibration). However, our algorithm has been designed to manage a huge quantity of very high-quality data (preselected), unlike SWAS [START_REF] Gao | SWAS: a shear-wave analysis system for semi-automatic measurement of shear-wave splitting above small earthquakes[END_REF], which allows processing of all the small earthquakes. As suggested by [START_REF] Crampin | A review of techniques for measuring shear-wave splitting above small earthquakes[END_REF], rejecting many records may lead to bias. However, due to the complexity of the waveforms collected in volcanic areas that are generally contaminated by the superposition of signals originated by the volcanic and/or hydrothermal dynamic, we preferred to apply severe selection rules on the analyzed dataset to reduce the risks of misinterpretations. Consequently, we measured the splitting parameters on signals showing clear shear wave onsets with high signal-to-noise ratios (S/N > 6) and incidence angles strictly inside the shear wave window (theoretically 35 • ). Moreover, all the selected data, after having been processed with the quasiautomatic analyses, passed the final hodograms check, so we did not reject any data at the final stage. As well known, an irregular surface topography severely distort the shear-wave window; we did not take into account this effect in our study, since all the stations we analyzed are installed in locally flat areas .

We found that at Mt. Vesuvius, the fast polarization direction strikes generally NS at BKN station and EW at BKE station; however, during the main seismic crisis which occurred in the last years (1996 and 1999), ϕ observations at the station nearest the seismic source (BKE) are rotated toward the NS direction (Fig. 3b). This result confirms previous observations by Bianco et al. (1998a[START_REF] Bianco | The seismic crisis at Mt. Vesuvius during 1995 and 1996[END_REF] and Del [START_REF] Pezzo | Changes in coda decay rate and shear-wave splitting parameters associated with seismic swarms at Mt. Vesuvius, Italy[END_REF] who analyzed the same dataset with different methodologies. The anisotropy percentage ξ is 4% for the background seismicity. Conversely, the same calculation performed only on data collected during the March-April 1996 and October 1999 swarms yields ξ = 8%.

In the Campi Flegrei area, the polarization of the fast wave shows an irregular pattern: for data collected during the 1982-1984 crisis, the qS 1 polarization is approximately NNW/SSE-oriented at W11 and W12 stations and NNE/SSW-oriented at W21 station; for the 2000 crisis dataset, the fast polarization at station SLF strikes NNW-SSE; for the 2004-2006 crisis the ϕ direction is approximately NNW/SSE-oriented at ASB2 and NE-SW at AMSE and BGNB stations. The spatial distribution of the qS 1 polarization appears to be complex (Fig. 4) and not easily interpretable. The anisotropy percentages span from 7% for dataset 1, to 4% for dataset 2, and 3% for the last dataset. The shear wave splitting parameters at Campi Flegrei have been investigated by [START_REF] Savage | Shearwave anisotropy of active tectonic regions via automated S-wave polarization analysis[END_REF]) and by [START_REF] Bianco | The role of hydrothermal fluids in triggering the July-August 2000 seismic swarm at Campi Flegrei (Italy): evidences from seismological and mesostructural data[END_REF]; both studies found results in agreements with the present paper.

At Mt. Etna, the fast polarization for all the investigated datasets strikes, on average, in the NS direction, which is the direction of the maximum compressive stress acting at regional scale [START_REF] Musumeci | Seismological evidence of stress induced by dike injection during the 2001 Mt. Etna eruption[END_REF]. ξ percentage is between 1% and 6%, the lowest value pertaining to the 1988 dataset collected in a period with no eruption, the highest obtained for data sampling the opening of an eruptive fracture during the 1989 eruption. Our results are compatible with those obtained by [START_REF] Bianco | Shear-wave polarization alignment on the eastern flank of Mt. Etna volcano (Sicily, Italy)[END_REF]Bianco et al. ( , 1998b[START_REF] Bianco | Shear wave splitting changes associated with the 2001 volcanic eruption on Mount Etna[END_REF].

Different circumstances may produce the measured shear wave splitting. We reject the hypothesis that S-P conversion was responsible for the observed splitting because at vertical incidence, this type of conversion shows the greatest amplitude on the vertical component [START_REF] Liu | Effects of the internal shear wave window: comparison with anisotropy induced splitting[END_REF], whereas the whole analyzed data exhibit maximum amplitude in correspondence with the splitting on the horizontal components. We exclude the influence of aligned mineral grains on our measurements since the investigated volcanic areas are not underlain by tilted sedimentary or metamorphic rocks including these aligned minerals, at least at the depth sampled by our raypaths [START_REF] Renzo | Magmatic history of Somma-Vesuvius on the basis of new geochemical and isotopic data from a deep borehole (Camaldoli della Torre)[END_REF][START_REF] Bousquet | The tectonics and geodynamics of Mt. Etna: synthesis and interpretation of geological and geophysical data[END_REF]. Waveforms that look like shear wave splitting (but are not) might be caused by complex source mechanisms [START_REF] Zollo | S-wave polarization inversion of the 15 October 1979, 23:19 imperial valley aftershock: evidence for anisotropy and a simple source mechanism[END_REF], but due to the low magnitude and small source dimension (presumed simple fault plane) of the events (M MAX = 3.6), we reject this hypothesis. For the above reasons, we suggest that the observed shear wave splitting is consistent with propagation through distribution of nearly vertical nearly parallel fluid saturated cracks striking in the direction of the neighboring system fractures. The same conclusions have been reached by the different authors that previously investigated most of the datasets analyzed in the present paper, as discussed in the other sections.

As suggested by [START_REF] Zatsepin | Modeling the compliance of crustal rock-I. Response of shear-wave split ting to differential stress[END_REF], for ξ ≥ 5%, a fracture criticality regime occurs. The measured values of anisotropy percentage may suggest that at Mt. Vesuvius, the fracture criticality regime has been reached during the 1996 and 1999 crisis; however, no eruption has been observed in these cases, the last eruption occurring on March 1944. At Campi Flegrei, the fracture criticality regime has been reached during the 1982-1984 crisis, and in this case too, no eruption occurred. Interestingly, the percentage of anisotropy in this area decreases with time; possibly, this observation may be explained with a lower degree of fracturing in the upper crust with a different circulation of hydrothermal fluids for the three crisis. At Mt. Etna, the fracture criticality regime was reached during the opening of the 1989 eruptive fracture, but not during the 2001 eruption, even though ξ assumes a borderline value. The correlation between reaching a fracture criticality regime and the occurrence of eruptions seems to be not straightforward for our datasets. The percentage of shear wave anisotropy and fracture criticality depend on several different factors such as Poisson's ratio, pore fluid velocity, and many others, being higher in areas of high heat flow [START_REF] Crampin | A review of the effects of crack geometry on wave propagation through aligned cracks[END_REF]. Our observations, on average, confirm this behavior, since we found normalized time delay generally greater than 10 ms/km. However, a general increase of ξ has been observed every time we analyzed data collected during volcanic crises in which significant movements of fluids have been documented (in particular at Mt. Vesuvius and Campi Flegrei, e.g., [START_REF] Bianco | The role of hydrothermal fluids in triggering the July-August 2000 seismic swarm at Campi Flegrei (Italy): evidences from seismological and mesostructural data[END_REF]Del Pezzo et al. 2004 and references therein). Our observations on the temporal variation of anisotropy percentages may account for a different regime in the fluid circulation favored by the critical opening of the stress aligned cracks, in agreement with the APE modeling [START_REF] Zatsepin | Modeling the compliance of crustal rock-I. Response of shear-wave split ting to differential stress[END_REF]).

Fig. 1

 1 Fig. 1 Examples of the two window lengths as determined by the program for Mt. Etna (a) and Campi Flegrei (b) events. The straight line marks the beginning of the window (corresponding to S-wave arrival time minus two samples); the dotted one indicates the end. The upper panel shows the window chosen by the program to obtain ϕ through

Fig. 2

 2 Fig. 2 The seismic network at Mt. Vesuvius; the location site of BKE and BKN stations are circled

Fig. 4

 4 Fig. 4 Rose diagrams of the fast polarization direction at Campi Flegrei area; stations and results for datasets 1, 2, and 3 are identified by the colors red, yellow, and green, respectively

Fig. 5

 5 Fig. 5 Map of the seismic stations that collected data on Mt. Etna analyzed in this study. The stations operated in different periods: DMT and POM operated in 1988 when no eruption occurred; NOC, CCV, and B92 collected data including seismicity that occurred during 1989 eruption; MNT and ESP operated during the 2001 eruption

Fig. 6 Fig. 7

 67 Fig. 6 Rose diagrams at Mt. Etna for the 1988 dataset at stations DMT and POM

Fig. 8

 8 Fig. 8 Rose diagrams at Mt. Etna for the 2001 dataset at stations MNT and ESP

Acknowledgements Thoughtful comments from David