Intrinsic Qp at Mt . Etna from the inversion of rise times of 2002 microearthquake sequence

About three-hundred microearthquakes, preceeding and accompanying the 2002-2003 Mt. Etna flank eruption, were considered in this study. On the high-quality velocity seismograms, measurements of the first half cycle of the wave, the so-called rise time τ, were carried out. By using the rise time method, these data were inverted to infer an estimate of the intrinsic quality factor Qp of P waves and of the source rise time τ0 of the events, which represents an estimate of the duration of the rupture process. Two kind of inversions were carried out. In the first inversion τ0 was derived from the magnitude duration of the events, assuming a constant stress drop and Qp was inferred from the inversion of reduced rise times τ−τ0. In the second inversion both τ0 and Qp were inferred from the inversion of rise times. To determine the model parameters that realize the compromise between model simplicity and quality of the fit, the corrected Akaike information criterion was used. After this analysis we obtained Qp=57±42. The correlation among the inferred τ0 and Qp, which is caused by some events which concomitantly have high τ0 (>30 ms) and high Qp (>100) indicates that the technique used is able to model rise time versus travel time trend only for source dimensions less than about 80 m. Mailing address: Dr. Salvatore de Lorenzo, Dipartimento di Geologia e Geofisica, Università degli Studi di Bari, Campus Universitario, Palazzo di Scienze della Terra, Via E. Orabona, 4, 70125 Bari, Italy: e-mail: delorenzo@geo.uniba.it


Introduction
The intrinsic quality factor Q p of the compressional body waves is considered one of the geophysical parameters best correlated to the physical state of the rocks.This is because, as shown in laboratory studies (Kampfmann and Berckemer, 1985;Sato and Sacks, 1989) a type-Arrhenius exponential law relates Qp to the temperature T and the pressure P of the rocks.In volcanic areas, low Q p values are usually associated with high temperature rocks (e.g., Sanders et al., 1995;de Lorenzo et al., 2001).However there are several reasons which make it difficult to determine the temperature of the rocks from estimates of Qp in the crust.First of all, Qp depends not only on the temperature but also on the percentage of fluid content in the rocks (Bourbiè et al., 1987).Morevoer, the Qp estimates are also dependent on the technique used to retrieve them (Tonn, 1989).These considerations lead to the conclusion that only a comparison between Qp and T in deep boreholes can help us to locally calibrate a Qp-T relationship (e.g. de Lorenzo et al., 2001) and then to infer the temperature field from threedimensional images.Some previous studies have examined the attenuation properties at Mt. Etna volcano by applying frequency domain techniques (Patanè et al., 1994(Patanè et al., , 1997)).The results obtained by the above studies have shown that the effect of attenuation on seismic radiation at Mt. Etna is considerable.Patanè et al. (1994) reported significant variations of Q from P waves as a function of depth.In particular, a drop of Q values for earthquakes located at depths less than 5 km was observed, supporting the idea that the upper part of the crust and the shallow volcanic layers are characterized by low Q.Recently a P-wave attenuation tomography study was carried out at Mt. Etna, down to 15 km of depth, using high quality data recorded in the period 1994-2001 (De Gori et al., 2005).Results from this first tomography were refined in the shallow layers (around 2 km) of the volcanic edifice by Martìnez-Arevalo et al. (2005), further confirming significant variations of total Q as a function of depth and azimuth.
In the frame of the European project named VOLUME (an acronym of «VOLcanoes Understanding of Mass movEment») aimed to investigate the movement of fluid masses inside Mt.Etna, a research line has been dedicated to the reconstruction of Qp images of the volcano with different techniques.This is because several studies have pointed out that the estimates of Q could depend on both the data to be inverted and on the technique used.The attenation parameter Q −1 estimated by the inversion of seismic spectra at Mt. Etna by De Gori et al. (2005) and Martìnez-Arevalo et al. (2005) is the sum of two parameters, Qp −1 and Qc −1 , which take into account two different and concomitant physical phenomena.Qp −1 is the intrinsic attenuation of P waves, which is related to the losses of energy due to the anelasticity of the Earth.Qc −1 is the scattering attenuation of the wave field and is related to the energy content of the secondary wave field, i.e. the wave field generated by the elastic heterogeneities inside the medium travelled by the waves (e.g., Mitchell, 1995;Del Pezzo et al., 2001).It has also been shown (e.g., Liu et al., 1994) that the methods based on the inversion of pulse widths of first P waves are more appropriate to estimate the intrinsic Qp, in that only the start of the signal, which is rela-tively free of complicated effects, such as complexities at the source and/or scattering from heterogeneities, is used.
Based on these reasonings, this article presents the results of a first study aimed to obtain a first estimate of the average intrinsic Qp at Mt. Etna by using the classic rise time method (Gladwin and Stacey, 1974).The analysis concerns about three-hundred microearthquakes recorded at Mt. Etna preceeding and accompanying the 2002-2003 Mt.Etna flank eruption.

The Mt. Etna geological framework
Mt. Etna is located in the complex geodynamic framework of Eastern Sicily, where major regional structural lineaments play a key role in the dynamic processes of the volcano (e.g., Patanè and Giampiccolo, 2003).
Several years of structural and geophysical observations have revealed that the orientation of most of the eruptive systems coincides with two structural trends NNW-SSE and NE-SW, observed both in the volcanic area and in the regional context (e.g., Azzaro and Neri, 1992).These alignments are hypothesized as the main volcanogenetic structures (e.g., Gresta et al., 1998) which control the evolution of Mt.Etna, as their interference establishes a weakness volume along which magma can rise from depth (Rasà et al., 1995).
Over the last 30 years Mt.Etna has had a high rate of eruptive events and therefore it constitutes one of the most important natural laboratories for the understanding of eruptive processes and lava uprising in basaltic-type volcanic environments.Seismic observations allowed us to carry out detailed investigations on major aspects of seismicity.In particular, the most recent effusive eruptions which occurred in 1989, 1991-1993, 1999, 2001 and 2002-2003 have offered good examples for quantitative analysis based on seismic data in digital format.
We focus on the 2002-2003 flank eruption which started in the night between October 26 and 27, 2002 with a seismic swarm in the central and upper part of the volcano.Fissures on both the NE and S flanks were activated with a huge lava emission and powerful explosive activity from the southern fracture field.The eruption ended on January 28, 2003 after 94 days.More than 800 events were recorded during this period by the permanent seismic network managed by INGV.Most of these events have magnitude less than 3.0 and only few earthquakes reached a magnitude duration M d = 4.4.

Technique
The rise time method is based on early experimental (Gladwin and Stacey, 1974) and theoretical (Kjartannson, 1979) studies.More recently, Wu and Lees (1996) showed that the theoretical relationship on which the method is based remains valid if we consider a heterogeneous anelastic structure.The most limiting assumption of this method is that it neglects the directivity effect of the seismic radiation generated by a finite dimension seismic source (Zollo and de Lorenzo, 2001).
On a velocity seismogram the rise time can be defined as the time interval between the onset of the phase and its first zero crossing time (fig.1).If the spatio-temporal finiteness of the seismic source process can be neglected, the variation of the rise time versus the travel time of the wave will be given by (3.1) where τ0 is the rise time at the source, Qp is the quality factor and T is the travel time.The constant C was found to be equal to 0.5 for a constant Q attenuation operator (Kjartansson, 1979).
The methods based on the inversion of the rise times are expected to give the most reliable estimates of the intrinsic attenuation (Liu et al., 1994).In fact, since only a very limited portion of the seismogram is used, the effects of multiple waves generated in the thin layers around the recording site are usually minimized.
Finally, as a consequence of the point-like source assumption, the Qp inferred by using eq.(3.1) can be considered as the minimum Q p estimated from pulse width inversion (de Lorenzo et al., 2004).

Data analysis
We used a data set of about 300 well located events (Erx, Ery and Erz <1 km, Gap < 100°a nd RMS < 0.20 s), with magnitude duration ranging between 1.4 and 4.4 and focal depth mainly concentrated in the first 5 km, recorded during an episode of intense seismic activity which preceded and accompanied the 2002-2003 Mt.Etna flank eruption.The earthquakes were selected from the catalogue among those re-located by using the most recent 3D velocity model (Cocina et al., 2005).Data from stations of the Mt.Etna permanent seismic network run by Istituto Nazionale di Geofisica e Vulcanologia of Catania (INGV-CT) were analysed.Most of the stations considered in this analysis are equipped with one-component short period (1 s) sensors, except 8 equipped with three-component sensors.Although other stations were operating during the occurrence of the analysed seismic sequence, data were not as well recorded as required for accurate seismic attenuation analysis in time domain.
We discarded the events for which the rise times are available at fewer than four stations.The remaining dataset was then composed of 147 events, all having typical tectonic-earth-   The origin time, the identification number and magnitude duration of the events used in the study, computed by INGV, are given in table I.A total of 1053 data was available for the study.
The rise time was measured on each first arrival P-wave by computing the time interval between the onset of the phase and its first zero crossing time.We discarded all data for which the P signal to noise ratio was lower than 10, and considered only the waveforms on which the onset of P-wave was clearly readable.
We did not perform the deconvolution for the instrumental response for the following three reasons: 1) First of all, the technique we are using does not need amplitude information.2) Second, the effect of filtering operated in the deconvolution for the instrumental response could cause the generation of unwanted artificial signals (see Mulargia and Geller, 2003 and references therein) which could bias the rise time estimate.3) Finally, if we exclude one datum, the observed rise times vary between 0.025 s and 0.3 s (fig.3), to which correspond an average frequency content ranging from a minimum of about 3 Hz to a maximum of about 40 Hz.Consequently the frequency content of the analyzed signals is always contained in the frequency band (1-50 Hz) where the instrument response is flat and cannot distort the duration and shape of the observed signals.
Only data not affected by multipathing during the first half-cycle of the wave were considered in the analysis.The multipathing effects  can often be easily recognized as a sharp discontinuity (e.g., de Lorenzo and Zollo, 2003) which breaks the approximately bipolar shape of the wave on the velocity seismogram (fig.4) and are generally caused by the presence of thin layers below the recording site, where part of the energy remains trapped and is subjected to multiple reflection.The quality of the onset varies with varying the level of noise of the recorded traces.For this reason we divided the dataset into three categories.A maximum error on rise time equal to 5% was estimated for data which have the highest quality.For data of intermediate quality a maximum error of 10% was estimated.Finally, a maximum error of 20% was associated to data having the worst quality.As a consequence an average error on rise time equal to 10 ms was estimated.
The plot of the seismic rays under the assumption of a homogeneous velocity model is shown in fig.5a,b.We infer that the centraleastern part of the array is very well illuminated by the seismic rays, until to a depth of about 3-4 km.The Qp estimates we present in the next section will be then particularly representative of the average Q p of this volume.

Results
We carried out two different kinds of data inversions.In the first inversion we assumed that the stress drop of the studied earthquakes is a constant, so that a self-similar behaviour of earthquakes with a different energy content is expected.Under this assumption the source rise time τ0 of each event was computed by using the seismic moment versus magnitude duration relationship calibrated for the Etnean area (Patanè et al., 1993) and the laws which relate source radius, seismic moment, stress drop, rupture velocity and source rise time for a circular crack model.In the second inversion, source rise times were directly estimated from data using the relationship (3.1).

Q p estimates under the assumption of a constant stress drop
For each of the 147 considered events, seismic moment M 0 was inferred from the magnitude duration Md (table I), using the relationship which relates M d to M0 for the Etnean area (Patanè et al., 1993) (5.1) with (5.2) . (5.3) For a circular crack, the source rise time τ0, which represents the time duration of the slipping on the fault, is related to the source radius L and the rupture velocity V r of a circular crack by the equation (Brune, 1970) .
(5.4) By combining eq. ( 5.4) with the relationship which combines stress drop ∆σ, M0 and L (Keilis-Borok, 1959) (5.5) we obtain . (5.6) In order to compute τ0, we need to estimate the rupture velocity V r of the earthquakes.Unfortunately, this parameter is poorly known for small magnitude earthquakes, owing to its correlation with the other source parameters (e.g., Deichmann, 1997).Theoretical and laboratory studies (Madariaga, 1976) indicate that Vr ranges between 0.6 V s and 0.9 Vs.Starting from the observation of a great dispersal of rise time versus travel time for the considered events, which is often attributed to significant directivity effects (e.g., de Lorenzo et al., 2004), we have assumed a high value of the rupture velocity (V r= 0.9 Vs=2.2 km/s).
For the Mt.Etna ∆σ is known to range from a few MPa to about 100 MPa (e.g., Patanè et al., 1994Patanè et al., , 1997;;Patanè and Giampiccolo, 2003).For this reason, first of all we carried out five different inversions by using five possible values of ∆σ (∆σ =1, 5, 10, 50, 100 MPa), and then we compared the results of the inversions.
To this aim, for each value of ∆σ, source rise times τ 0 of the events were computed using eq.(5.6).Then, for each event, the reduced rise times τ−τ 0 were inferred from eq. (3.1) and inverted to estimate Qp.
At the first step of the inversion we assumed a constant Q for the entire area.The inversion of the entire data set of rise times was performed using a weighted linear inversion scheme with weigths equal to the inverse of the variance of data.The results are summarized in fig.6a-c.It is worth noting that about the same quality of the fit is obtained for ∆σ =10 MPa (an average residual of 42 ms in L1 norm and an average residual of 62 ms in L2 norm) and for ∆σ =50 MPa (an average residual of 42 ms in L1 norm and an average residual of 61 ms in L2 norm).The coefficient of correlation for the two cases, as computed using the relationship for weighted data (Green and Margerison, 1978) is also about the same (R 2 = 0.24).For this reason, Qp was estimated by averaging the two Qp estimates for ∆σ =10 MPa (Qp∼33) and ∆σ =50 MPa (Qp∼30), obtaining . (5.7) Owing to the difficulties in accurately accounting for both the errors on data and the tradeoff among source parameters, instead of fixing the stress drop to one of the values which produce the comparable fit to data (∆σ =10 or ∆σ =50 MPa), we preferred to consider, in the calculation of Q p of each event, both the cases ∆σ =10 MPa and ∆σ =50 MPa, choosing, for each event, the Qp value which gives rise to the best fit result.Results are summarized in table II.Only two events show a negative non physical Q p value: event #85, for which we have Qp=−37 and event #100, for which we have Qp=−100.Moreover, only event #21 shows a very high value (Q p=1230).For the remaining events the Qp variations are more limited, with a minimum Qp=9 and a maximum Qp=96.The weighted average of Q p, using the results of all the 147 inversions and as weights the inverse of the residuals in L1 norm is exactly the same as the Q p obtained at the first inversion step . (5.8)

Joint estimation of τ0 and Qp
In another attempt we inverted the rise times to jointly infer τ0 and Qp of each event.
In a first inversion run, we used a linear weighted inversion scheme, with weights equal to the inverse of the variance of data.Unfortunately this approach does not allow us to impose positivity constraints on τ 0 and Qp.For this reason, after the inversions, only for some events we inferred positive values of both τ 0 and Qp, reported in table III.To overcome the problem of negative values of τ0 and/or Qp for the remaining events, we used a non linear in- Table III.τ0 and Qp estimates of the events; the lines with bold characters indicates the events for which the Simplex method has been used to overcome the problem of negative τ0 and Qp values.version scheme based on the use of the Simplex Downhill method (Press et al., 1989) which allowed us to impose positivity constraints on both τ0 and Qp.To estimate the error on the inferred model parameters we applied a statistical approach based on the use of the random devi-ates (Vasco and Johnson, 1998).This consists of computing N rand datasets by adding to each data a random quantity selected in the range of the error affecting them.The inversion results of the N rand datasets can then be used to estimate the average values of the model parameters and (5.9) A further variance reduction was obtained with respect to the previous case of a heterogeneous Qp but a constant stress drop ∆σ.In fact the average residual in L1 norm is now equal to 23 ms (a residual reduction of 33%) and the average residual in L2 norm is equal to 37 (a residual reduction of 30%).The trend of residuals versus travel times is shown in fig.7b.

Application of the «Occam razor principle»
One point that needs to be addressed concerns the evaluation of the statistical significance of the two models used to fit data.Since the two models are characterized by a different number of model 68 51 Qp != parameters it is important establish, for each event, what of them satisfy the so called Occam's razor principle, i.e. realizes the best compromise between the quality of the fit and the simplicity of the model.To solve this problem we used the corrected Akaike Information Criterion (AICc) (Cavanaugh, 1997;Cavanaugh and Shimway, 1998), which represents a modified version of the Akaike (1974) criterion, to account for incompleteness of data.This consists of comparing, for the two models, the value of the parameter (5.10)where k is the number of model parameters, n the number of data and L represents the likelihood function, given by .(5.11)The model which gives rise to the minimum AICc has then to be considered the most significant from a statistical point of view.
Table IV summarizes the comparison between the AICc values for the two considered models.In this way we were able to infer that, for 59 events the above choice of retrieving source rise times assuming a constant stress drop produces a better fit of data, whereas for the remaining 68 events the source rise times are better constrained from the inversion of rise times.By averaging the Qp estimates obtained with this analysis we finally obtained . (5.12) Figure 8 shows the trend of the Qp of each event estimated with the AICc as a function of τ0.It is worth noting that a residual correlation among τ0 and Q p is inferred.42 57 Qp !=

Discussion and conclusions
As thoroughly debated in Wu and Lees (1996) and de Lorenzo (1998), the slope C of the straight line interpolating τ versus Q p -1 depends on the frequency content of the source; by calibrating eq.(3.1) with different source time functions, Wu and Lees (1996) showed that C=0.5 has to be adopted if it is assumed that the the signal generated at the source has a Gaussian shape, which is a smooth representation of the unipolar source time function generated by a circular crack.By taking into account that a more realistic source cannot have the same low-frequency content as the Gaussian function, Wu and Lees (1996) showed that the C could be slightly higher.They estimated as  1230 Let us consider now the residual correlation between Qp and τ 0 .First, we note that the Q p es-timate obtained using the AICc (Q p= 57±42) is intermediate between the Qp estimate obtained under the assumption of a constant stress drop and that obtained by jointly retrieving τ0 and Qp.However, this result does not necessarily im-ply a departure from the constancy of the energy release per seismic moment, nor is it necessarily caused by a great heterogeneity in the attenuating properties of the area, as one can suppose by taking into account the high value of the standard deviation on Qp.A possible interpretation of the correlation between Q p and τ0 may be that, in many cases, the data available are not sufficient to correctly estimate the two parameters.As an example, it can be easily demonstrat-Fig.10.Plot of rise time versus travel time for some studied events.On each plot the best fitting straight line, obtained from the inversion of rise times, is superimposed.ed that a strong clustering of empirical data in a small range of values for both the considered variables induces a correlation between the regression parameters.Thus the uneven and insufficient sampling for many events could be responsible for the observed phenomenon.
Alternatively, we have to account for the limitations of the hypothesis on which the used technique is based.In fact the adopted technique uses a very simplified source model, i.e. an isotropic seismic source.It is known that, for a finite dimension seismic source, the directivity source effect on the rise time of first P waves can be described through the relationship (Zollo and de Lorenzo, 2001) (6.1)where θ is the takeoff angle (the angle formed by the tangent to the ray leaving the source and the normal to the fault plane) and Vp the P-wave velocity of at the source.From this equation it immediately follows that the variations of the source rise time around the non-directive source rise time given in eq.(5.4) increase with increasing the source dimension L.Moreover, for a directive source time function, the equation which describes the relationship between τ and Q −1 is non linear, as demonstrated in Zollo andde Lorenzo (2001) andde Lorenzo et al. (2004).In particular Zollo and de Lorenzo (2001) showed that the higher the source dimension is, the higher will be the difference between the rise times predicted by the directive source model and those predicted by eq.(3.1).In our case, all the data which are responsible for the observed correlation among τ 0 and Qp (the data shown in the grey square of fig.8) have been selected by AICc among the model parameters obtained from the joint inversion of τ0 and Qp, i.e. without imposing a constant stress drop for the area (fig.9).This could imply that, above a given threshold (τ0~30-40 ms) the directivity source effect tends to be so high to produce a non linear trend of rise times versus the travel time, which may result in a bad estimation of the estimated τ0 and Qp if one uses the linear interpolating scheme given by eq.(3.1).This result is also supported by the visual inspection of the quality of the matching of the model to data, as shown in fig. 10 for some events.It is quite evident that, even if, in some cases, an increase in rise time versus travel time can be observed, the dispersal of data around the best fit straight line is significantly higher than the error on data.A possible bias in Q estimates could be caused by the different instrumental responses.It has to be noted that twenty-two of the available stations are equipped with Mark L4C, having an eigenfrequency of 1 Hz.The remaining two stations are equipped with Lennartz Le3D, having an eigenfrequency of 0.05 Hz.For both the two kind of stations the cutoff frequency is around 50 Hz.Since the average frequency content of the pulses ranges from 3 to 40 Hz, i.e. lies in the linear part of the amplitude and phase response of the instrument, we do not expect a significant distortion in amplitude of signals and phase shifting, even if the signals having the higher average frequency content (around 40 Hz) could be more sensitive to the non linear part of the instrumental response.However, only a dozen of data (fig.3) have an average frequency higher than 30 Hz and then this problem affects only a very negligible portion of the dataset and cannot affect the obtained average Q estimate.
Therefore we conclude that, for the Etnean area and in the case of a limited exploration of travel times, the applicability of the classic linear rise time method has to be limited only to those earthquakes which have a source rise time less or equal to about 35 ms, at which corresponds, using a rupture velocity Vr=0.9 Vs=2.2 km/s, an average source dimension of about 80 m.

Fig. 1 .
Fig. 1.Schematic picture of rise time τ as measured on a velocity seismogram.

Fig. 2 .
Fig. 2. Epicenters of the events (grey circles) and recording stations (black triangles) considered in this study.

Fig. 4 .
Fig. 4. A velocigram of a microearthquake occurred at Mt. Etna.The discontinuity occurring during the first half-cycle of the first P-wave, due to a secondary arrival, impedes the measurement of the rise time on this seismogram.

Figc
Fig. 6a-c.a) Results of the inversions of rise times for different values of ∆σ, assuming a constant ∆σ and a homogeneous Qp.On the top-left, the plot of the average residual in L2 norm versus ∆σ is shown.On the topright the plot of the average residual in L1 norm versus ∆σ is shown.On the bottom-left the plot of average Qp versus ∆σ is shown.On the bottom-right the plot of the squared correlation coefficient versus ∆σ is shown.b) Plot of reduced rise times versus travel times for ∆σ=10 MPa and ∆σ=50 MPa.c) Plot of rise time residuals versus travel times for ∆σ=10 MPa and for ∆σ=50 MPa under the assumption of a homogeneous Qp.
However, with respect to the previous result, based on the assumption of a homogeneous Qp, the assumption of a heterogenous Q structure gives rise to a further variance reduction.In fact the average residual is now 34 ms in L1 norm (a residual reduction of 19%) and 53 ms in L2 norm (a residual reduction of about 13%).The Q 32 2 p ! = trend of residuals versus travel time is shown in fig.7a.

Fig.
Fig. 7a,b.a) Plot of rise time residuals versus travel times under the assumption of a constant ∆σ and a heterogeneous Qp. b) Plot of rise time residuals versus travel times after the joint inversion of τ0 and Qp.

Fig. 9 .
Fig. 9. Plot of τ0 versus the identification number of the event and of Qp versus the identification number of the event.

Fig. 8 .
Fig. 8. Plot of Qp versus τ0 of the studied events after the application of the AICc.The points enclosed in the grey rectangle are mainly responsible for the observed correlation.

Table I .
The dataset available for this study.Md = magnitude duration; Id = identification number of the event; Date = month-day-hour and minute of occurrence of the event (year: 2002).
Intrinsic Qp at Mt. Etna from the inversion of rise times of 2002 microearthquake sequence quake waveforms, whose localizations, deduced from the INGV catalogue, are shown in fig.2, together with the position of the recording stations.

Table II .
Qp estimates under the assumption of a constant stress drop.