R Davi

G S O'brien

L De Barros

I Lokmer

C J Bean

Philippe Lesage

Mauricio M Mora

Gerardo J Soto

G S O"brien

Seismic source mechanisms of tremor recorded on Arenal volcano

Keywords: Arenal volcano, tremor, moment tensor inversion, source mechanisms Arenal volcano, tremor, moment tensor inversion, source mechanisms

HAL is

Introduction

Volcanoes are generally characterised by a wide range of seismic events. This variety is normally due to complexity within these types of environments where different phenomena such as fluid movement, rock fractures, heterogeneity, topography, and near surface velocity structures can generate and affect the seismic response. Volcano-tectonic (VT), long-period (LP) and tremor are some examples of seismic records typical of volcanoes (McNutt, 2005). While the source mechanisms generating VT events are generally well-known, LP and tremor source mechanisms are not fully understood. The generation of these types of events is currently thought to be associated with fluid movements (Chouet, 2003), but a full understanding of these source processes is still very difficult to obtain. Since these types of events may be used as precursory phenomena for volcanic eruptions, a full understanding of their generative mechanisms is of crucial importance to constrain the dynamics of the plumbing system of restless volcanoes. One of the most common tools to derive the seismic source mechanisms is a moment tensor inversion. This approach assumes that the moment tensor components can be represented as a system of equivalent forces that produces the same wavefield as the actual physical processes acting at the source. Moment tensor waveform inversions have been applied to a wide range of volcanoes around the world (Chouet et al., 2003;Nakano et al., 2003;Kumagai et al., 2005;Lokmer et al., 2007;De Barros et al., 2011). However, while inversions for very long period events (VLP) have been successfully performed (Ohminato et al., 1998;Chouet et al., 2003), the shorter wavelengths of LP events are generally distorted by strong topography and structural heterogeneities. These distortion effects can influence the inversion results, leading to erroneous solutions (Bean et al., 2008). Although the problems related to the inversion of LP and VLP events have been widely studied (Nakano et. al., 2007;Cesca and Dahm, 2007;Bean et al., 2008;Davi et al., 2010), currently there have been no attempts to invert tremor. This is probably due to the difficulty in isolating a clear first onset of the tremor waveforms, the importance of which derives from the idea that the initial part of the signal can be considered to represent a source effect only, while over time the signal becomes more influenced by path effects.

Furthermore, the long duration and the complexity of tremor make the inversion procedure, and consequently retrieving the source mechanism, a challenging task. In this article, we present the results of the inversion of 15 tremor events selected from a dataset recorded on Arenal volcano (Costa Rica) during a seismic experiment carried out in 2005. At Arenal, it is possible to isolate the first onset of tremor, hence these events provide a good opportunity to broaden our knowledge about tremor source mechanisms. The source position is determined by performing moment tensor grid searches within a volume of possible source points located under the crater summit. A simultaneous grid search for the source mechanisms and location is performed for each event. The source location is identified through the evaluation of the misfit values where the lowest value is considered the best estimation. Subsequently we defined a common source position evaluating the joint probability density function of the misfit values for all the events. For this common source point, solutions for moment tensor components (MT) and moment tensor components plus single forces (MT+SF) were evaluated. As shown later in the text, it is sensible to assume that all the tremor signals were generated at a common source location. The solutions have been analysed using a sliding time window over the retrieved source time function in order to closely analyse the recovered mechanisms. This procedure is assessed using synthetic tests in order to evaluate the reliability of the inversion procedure when the sliding time window is used across 5 minute long synthetic signals. The same procedure is then applied to the real data allowing us to investigate possible temporal variations or stability of the source mechanism generating the tremor events in Arenal.

Results show temporal stability of the source process generating the selected events. Furthermore, decomposition analysis is performed, which involved analysing the maximum amplitude of the signal considering the polarity of each moment tensor component; these results indicate the presence of a strong and constant isotropic component in the solutions.

Arenal and tremor models

Arenal is a 15 km 3 andesitic strato-volcano located in north western Costa Rica and rises to 1750 m a.s.l., 1100 m elevation with respect to its base (Alvarado and Soto, 2002, figure 1). Since 1968 Arenal is characterized by permanent activity accompanied by gas emissions, Strombolian-Vulcanian explosions, lava and pyroclastic flows. Arenal"s seismicity is mainly characterized by long period events, explosions and tremor. All these kinds of event, including tremor, can be accompanied, or not, by ash emission and sometimes ash can be emitted aseismically (Valade et al., 2011). Tremor events are the most common type of event in Arenal with durations up to several hours per day. Depending on the frequency content, two different types of tremor, spasmodic and harmonic, can be identified in Arenal, where spasmodic tremor is usually characterized by a wide range of frequencies from 1 to 6 Hz while harmonic tremor spectra contain evenly spaced peaks with fundamental frequency in the range [0.9-2] Hz (figure 2). They both display large temporal variations with gliding phenomena (positive and negative) and jumps in frequency that correspond to sudden decreases following positive gliding (Lesage et al., 2006;figure 3). Different models have been proposed to explain the origin of tremors in Arenal. These models can be divided into three main categories.

According to Benoit and McNutt (1997), tremor events are produced by a fluid-filled resonator (like an organ pipe), where the source fluid is a gas-charged magma. The organ pipe model could be used to explain the evenly spaced peaks of the harmonic tremor. Small variations in the gas content produce strong velocity changes that explain the observed frequency fluctuations. The excitation of resonance is not included in the model.

A second category includes models in which tremors are thought to be generated by the non-linear excitation of a visco-elastic channel caused by different pressure states of a fluid flowing inside the conduit (Julian, 1994;2000). The non-linear excitation is used to model changes in fluid-flow that could account for some of the seismic signals observed at Arenal. These models are based on the assumption that the high velocity fluid is flowing through conduit constrictions. However according to Lesage et al., (2006 and references therein) using a value of a few meters for the conduit diameter and considering the low effusive rate of about [0.3-0.6] m 3 s -1 the resulting magma flow velocity is of the order of a few millimeters per second. Therefore the physical conditions that could produce a non-linear excitation are not fulfilled in this volcano (Lesage et al., 2006).

One of the most recent models (Lesage et al., 2006) proposes that tremor events can be generated in a conduit filled with bubbly magma, where the upper part is blocked by solid or very viscous lava acting as a barrier for the gases. The plug is fractured in a way that allows gases to escape. The tremor is generated by repetitive pressure pulses that interact with standing waves in the conduit, a similar process as a clarinet. A detailed analysis of several mechanisms that can generate seismic waves by a fluid flow in a narrow channel suggests that the clarinet model is the most plausible one in volcanoes (Rust et al., 2008). The behaviour of the system is closely related to the properties of the lava plug. These properties are continuously changing with the intrusion of fresh magma into the lava pool and on the way the lava cools down and solidifies to create the plug. The later model is also consistent with joint observations of explosions and tremors by seismometers and Doppler radar (Valade et al., 2011). Their analysis shows large variability of the phenomenon and non-systematic relationship between seismic and radar signals. Tephra emissions were observed during part of the explosions and also during some episodes of tremor and seismic quiescence (Valade et al., 2011).

In this paper, we intend to investigate further and provide additional information on the seismic source mechanisms generating harmonic tremors, in order to broaden our knowledge about the plumbing system.

Data and experimental set-up

In February 2005, a seismic experiment was carried out on Arenal volcano by the University of Costa Rica and the Université de Savoie in France. During the experiment ten broadband seismic stations were deployed around the volcano. Three-components Guralp CMG40T 60-secondsseismometers with mini-Titan recorders were operating in continuous mode at 125 samples per second. Arenal is a dangerous environment due to the strong explosive activity, therefore it is not possible to access the top of the volcano since 2000. The hazardous environment, together with the steep topography and the rich vegetation on the eastern flank of the volcano, strongly influenced the station distribution adopted in the experiment (figure 1). The highest station, ARE1 (at ~750 m a.s.l) is located on the western flank at roughly 2 km from the summit. A pressure sensor and a microphone were also positioned in the same location. The experiment lasted 10 days (from the 11 th to the 21 st of February), during which a wide range of seismic events was recorded (e.g. explosions, LP events, and tremor). The main energy of the different events ranges between 0.8 and 5 Hz. The most common type is tremor, which can vary in duration up to several hours. Often the tremors appear to be triggered by an explosion quake, however the triggering mechanism is not always clear.

We selected 15 harmonic tremor signals with fundamental frequency of 0.8 to 1.5 Hz. The complex behaviour of the volcano makes the analysis of the dataset a difficult task. In addition, due to instrumentation problems, the number of recording stations is not constant for the whole length of the experiment, changing from a minimum of seven to a maximum of ten. For each of the selected tremor events, deconvolution of the instrument response is performed. Filtering, spectrum and spectrogram analysis have been carried out in order to gain a better understanding of the tremor events. Their spectrograms are characterized by spectral lines at different frequency intervals. The evenly spaced overtones are interpreted as resulting of a Dirac comb effect due to a stable repetitive source process (Lesage et al., 2006). Moreover, the spectra are made more complex by site and propagation effects due to heterogeneities in the medium and strong topography (Mora et al., 2001;2006). For example, while the frequency peaks are the same at different stations, strong amplitude variations of the spectral peaks can be observed between close stations. For each of the 15 tremor events, the whole length has been selected. The duration of the tremor is not constant for all the events and it varies from around 100 seconds to more than 500 seconds. Mean and trend have been removed from the dataset. A cosine taper and a bandpass filter (0.8-1.5 Hz) were applied prior to the inversion in order to isolate the fundamental frequency.

Location of tremor events

The location of the seismic source is crucial to gain a better understanding of the volcanic system and to improve the eruption forecasting (Métaxian et al., 2002). Tremor events are generally characterized by emergent onsets with waveforms affected by scattering effects due to the propagation through a highly heterogeneous medium and the interaction with the topography. All these characteristics make the use of the classical location methods (such as particle motions or similarities between waveforms) ineffective and imprecise. For this reason, in this paper, the location of the source is estimated through the probabilistic approach (described further in the text) using grid searches performed for the whole data length of 15 tremor events. We created a volume of 4735 possible source points (with a grid step of 15 m), located under the crater summit (480 x 380 x 840 m 3 ), the dimensions of which have been retrieved in previous studies of source location (Benoit and McNutt, 1997;Hagerty et al., 2000;Métaxian et al., 2002;2009). For each tremor event, a moment tensor grid search is performed in order to retrieve the waveform misfit values (R) for each point of the source volume and for the solutions with and without single forces. For the linear inverse problem d=Gm, the misfit between observed (d) and calculated data (Gm) in the least square sense can be expressed by:

    d d Gm d Gm d R T T    (1)
where d , G , m are the data, the Green"s functions and their derivatives, and the moment tensor and single force components, respectively. The superscript " T " indicates the transposed matrix.

The lowest misfit is considered the best estimation.

The single forces can be either considered (the solutions are moment tensor and single forces, MT+SF) or omitted (solution are only moment tensor components, MT). The Green"s functions and their derivative were calculated over the above mentioned source zone of 4735 points using 3D full wavefield numerical simulations (O"Brien and Bean, 2004). The topography and 1D velocity model are included in the simulations. The one-dimensional velocity was retrieved from sounding using the spatial autocorrelation method (SPAC) and seismic refraction experiments carried out on Arenal in 1997 (Mora et al., 2006). The numerical domain consists of a 13 x 11 x 6 km 3 where the topography is derived from the Digital Elevation Model (DEM) of the volcano using a spatial grid step of 15 m. The Green"s functions and their derivative are calculated using the Reciprocity theorem (e.g. Aki and Richards, 2002). The method and the inversion techniques are described in Davi et al. (2010). The result of this calculation is a set of 15 second long Green's functions with a frequency of up 5 Hz.

We analysed the time delays between the signals recorded at different stations and found that these delays are equal for each event. These results suggest a common source location for all the events.

With this assumption we performed grid searches for each tremor event, results of which indicate a point of minimum misfit for both solutions (moment tensor, MT and moment tensor plus single forces, MT+SF, respectively). To find a common source point for all the tremor events, we evaluate the joint probability density function P(dj) i.e. the likelihood that the common source is located at a certain position of the considered volume (Menke, 1984):

           N = i j i j ) (d R d P 1 2 1 exp (2)
where d j denotes the j-th source point, and N is the number of events. In order to find a common source position, inversions for MT and MT+SF for each of the 15 tremor events have been calculated, and eq. ( 2) is then applied to find the joint probability density function P(d j ). By evaluating the maximum of the joint probability density function, we were able to locate the common source position in an area of the source volume located under the crater summit for MT solution (the minimum misfit value for MT+SF solution is located nearby). Due to the similarity of the MT and MT+SF solution, and since Davi et al. (2010) stated that in the case of Arenal and with this station configuration the MT solution is more stable, we decided to fix the source location in the position retrieved by the grid search for the MT solution (figure 4, panel A). Figure 4,panel B shows the area of the volume with a 90% percentage of probability of finding the real location of the source. Since this area is not large (it has an extension of 100 m in each direction) we can confirm the reliability of our source location. The source is located 100 meters below the crater summit. The results of the grid search indicate that the retrieved position is in agreement with the results of Métaxian et al. (2002;2009). According to these authors, the source position of the tremor is located in a small area centred under the crater summit.

Synthetic tests

The inversion of tremor can be a sensitive procedure; in fact the complexity and the long duration of these types of events, along with errors in the velocity model, can easily influence the inversion results. For these reasons particular care must be taken in the interpretation of the results. Davi et al. (2010) performed some synthetic tests in order to verify the consistency and limitations related to the inversion procedure on Arenal volcano. The tests were performed using the best estimate of the velocity model and a station configuration which mapped the real location of ten stations deployed around Arenal volcano in 2005. The results highlighted how spurious single forces are easily generated under conditions common on volcanoes, such as noisy data, imperfect knowledge of the velocity model, and mislocated source positions. These phenomena, common of volcanic environments, can influence the correct retrieval of the source time function and lead to erroneous mechanisms. This information was used in the interpretation of the results of an inversion for an explosive event recorded on Arenal in 2005. Since tremor events can be seen as the product of a series of repetitive source time functions and, thus, differing only in their longer duration, the same results are considered valid for the inversion of the tremor signals. Hence, only the moment tensor components for moment tensor (MT) and moment tensor plus single forces (MT+SF) solutions were analysed while the single forces were neglected (see also De Barros et al., 2011).

Compared to the previous study (Davi et al, 2010), an additional difficulty in inverting tremor is the length of the signal. It is computationally very expensive to invert signals with duration of hours.

Hence by using synthetic data we can derive the best strategy to find the source mechanism of long lasting signals. For these reasons, we tested the inversion procedure designing two different types of synthetic tests in which (i) we inverted for the whole length of a synthetic tremor in order to test how the inversion procedure performs when long windows of data are utilized, and (ii) performing the inversion using only the initial part of the signal. The solutions were analyzed using a short sliding window (of 32 s) in order to verify the stability or variability over time of the mechanism generating the tremor events. This analysis can be made using the maximum amplitude of the signal (denoted MAS in this article), taking into account the polarity of each component (Chouet et al., 2005). An additional test is performed to investigate the sensitivity of the inversion to a poorly detailed velocity model.

From the convolution between a Dirac comb and a Ricker wavelet function with a central frequency of 2Hz, we created a 5 minute long synthetic tremor signal with the same spectral characteristics as the recorded signals. The source position is in the location point found through the grid searches.

We simulated two source mechanisms (i) a pure volumetric (M xx =M yy =M zz =1; M xy =M xz =M yz =0) and (ii) a compensate linear vector dipole (CLVD) mechanism (M xx =M yy =1; M zz =-2, M xy =M xz =M yz =0).

Inversion of the full length of the signal

In the first test an explosive mechanism was simulated using a 5 minute signal. The longer the duration of the tremor, the larger the required computational time. However the solutions for the entire length appear to be very stable with a misfit value approaching zero (figure 5). With a perfect velocity model the source time function and mechanism are perfectly retrieved. We selected a 32 second long time window and we slide this window through the solution. For each time window we decompose the solutions of MT and MT+SF components. The decomposition results of the pure volumetric mechanism show a [1 1 1] ratio of the moment tensor diagonal components (M xx , M yy , M zz ) in every interval for both solutions, MT and MT+SF. The percentage of isotropic components (% ISO) is roughly 100%, while the double couple (DC) and compensate linear vector dipole (CLVD) percentages are very close to zero. Results are shown in figure 6. The same test (not shown here) was performed using the pure CLVD mechanism [1 1 -2] and again the results of the decomposition indicate a 100% of CLVD and a 0% of ISO and DC components. When perfect conditions (correct velocity model and the whole length of signal) are used, the mechanism and the source time function are perfectly retrieved.

A second test was performed in order to assess the role played by poorly detailed velocity model when long and complex waveforms are used in the inversion procedure. Due to lateral heterogeneities and potential scatters the real velocity model is commonly more complex than what we use for modelling. Consequently, to resemble a real-like situation in order to test the influence of the velocity model on the inversion results we used a set of Green's functions calculated for a simpler velocity model (e.g. a homogeneous velocity model), to invert a signal calculated for a more realistic model (the layered model derived by Mora et al., 2006). We simulated a geometry of an horizontal crack and we performed the inversion using the 5 minutes synthetic signals. Results of the inversion for MT and MT+SF are shown in figure 7 (panel A). The solutions have been subsequently decomposed using again a window of 32 seconds. Figure 7 (panel B) shows the decomposition results. The solution without including the single forces (MT) appears to be more reliable showing a dominant isotropic component within ~65%, which remains constant along the solution. From these results we can conclude that, despite using Green's functions which are simpler than the true model, the relative percentage of ISO seems to be correctly retrieved. On the contrary the deviatoric part of the solution shows more inaccuracies. However the relative percentages of ISO, CLVD and DC remain constant along the solution which further validate the use of our decomposition method in order to assess the stability or variability of the source mechanism.

Inversion of partial length of the signal

. As the whole length of the signal cannot be always inverted because of the computation cost, the initial part of the signal is chosen because this is considered to be less affected by scattering and The observed inability to obtain a correct solution for the final part of the incomplete signal can be easily explained if the tremor signal is seen as a set of consecutive sources (like a Dirac comb).

Since each of these consecutive sources has a propagation effect of the same length as our Green"s functions, only those sources for which the whole length of the propagation effect is included in the signal can be successfully retrieved. In other words, the retrieved source-time function will be correct within a time window which starts at the beginning of the signal and whose length is equal to the difference between the length of the signal and the length of Green"s functions.

In summary, in ideal cases when the entire length of the signal is used, the correct source mechanism and source time function are retrieved for the whole length of the signal. This allows us to analyse such a solution within multiple sliding time windows in order to verify the stability or the variability of the mechanism generating the tremor. Hence the use of the whole length of the signal will allows us to detect, if present, variations in the source mechanisms and to analyse the influence of propagation effects through the signals. On the contrary, when inverting for a part of the signal, only the beginning part (length of the data minus length of the Green's functions) of the retrieved source-time function can be considered as a correct solution. Therefore this approach can be useful for inverting very long tremors, when the full length of the signal cannot be included in the inversion procedure.

Tremor inversion results and discussion

We selected 15 real tremor events recorded in Arenal volcano and performed a moment tensor inversion to retrieve the source mechanisms generating these events. The waveform inversion of the tremor is performed in the frequency domain for moment tensor solutions with and without single forces. For each of the selected events, the whole length of the tremor has been used to perform the inversion; the shortest events have durations of 100 s while the longest ones last for 600 s. All the selected tremor signals show a clear but emergent onset that can be easily isolated. From the analysis of the data, it is clearly noticeable that the waveforms differ from one event to another in length and amplitude. However the solutions of the inversions for each tremor show very similar patterns. An example of the moment tensor solution is shown in figure 9 in which the obtained MT+SF and MT solutions are plotted for the event recorded on February 13 th at 13:00. The fit between calculated and observed data is illustrated in figure 10. The fit (R= 0.6) shows a reasonable match for nearly all the stations which is indicative of the stability and accuracy of our inversion solution. In some cases the difference in amplitude between the observed and calculated data can be originated by local site effects (see station 4). These results are similar for the 15 events. However, in the MT+SF solution large single forces, especially in the z direction, are generated. From the results obtained by the synthetic tests (Davi et al., 2010), we highlighted how, in the case of Arenal, the errors caused by the poorly constrained velocity model, and the non-ideal station distribution, easily generated spurious single forces. Hence, we cannot easily distinguish, when inverting the real dataset, if, and to what extent, these components are real or spurious. For these reasons we excluded the single forces from the decomposition procedure. The same conclusions were reached by De Barros et al. (2011) for Mt. Etna volcano, Italy. We selected a constant window of 32 seconds and performed the decomposition sliding the window over the whole length of the solution for each event. For each window the percentage of isotropic, compensate linear vector dipole and double couple are calculated in order to retrieve the source mechanism generating the events. Results of the decomposition are shown in figure 11 for MT and figure 12 for MT+SF solutions (shown here for comparison with the MT solution). For each tremor event a stable isotropic percentage (ISO%) is found. The ISO% varies from 40% to 60% with some sporadic examples at 30%. Therefore the ISO% shows a small variation with a standard deviation in the order of the 5% from the mean value.

These results are common for the moment tensor components of both solutions, with and without including the single forces. The CLVD and DC components show a more scattered pattern, especially for the solution including the single forces. The CLVD component shows a standard deviation in the order of 12% while it reaches 17% in the case of the double couple components.

The same results are obtained (not shown here) when shorter and longer (e.g. 16 seconds and 40 seconds, respectively) decomposing time intervals are utilized. Results were also tested with an overlapping window (length of 24 s overlapping every 12 s). No significant variations were found.

The stability of the retrieved source mechanisms highlights the reliability of our solutions which do not seem to be strongly affected by path effects.

As shown by the synthetic tests, incorrect source mechanisms can be retrieved when there are significant uncertainties in the data. Since the level of uncertainty increases (due to poorly detailed velocity models, heterogeneities, network configuration, or source location errors) in the case of real waveform inversions, the percentage of these incorrect source mechanisms can increase. Since the decomposition procedure of the deviatoric part of moment tensor (i) is not uniquely defined and (ii) can be very sensitive to the noise contaminating MT solutions, the interpretation of the mechanism can be ambiguous. Furthermore, the outcome of the decomposition confirms that, in the case of Arenal volcano, with this station configuration and with the velocity model we employed, the solutions without the single forces are more stable than the MT+SF solutions. The same conclusion was reached previously (Davi et al, 2010) when the inversion procedure was tested for explosive events. For this reason we will base the interpretation of our results on the evidence shown by the MT solutions.

The decomposition method indicates solutions with 50 ± 5% of isotropic component, 35 ± 12% of CLVD and 15 ± 17% of DC for all the events. The overall trend shows a solution which is very stable in time and from one event to another. This confirms that the same non-destructive source processes generate the tremor and that the mechanisms do not vary in time. The diagonal components are, for each event, in a ratio of roughly [1 1 2]. Since the deviatoric part is not perfectly solved and our solutions show large uncertainties in the DC percentage, these results could be interpreted as the mechanism of a crack that is opening in a direction non-perpendicular to the crack plane. The analysis of the eigenvectors (figure 13) shows that a sub-horizontal crack with constant orientation is found for the 15 tremor events. The orientation of the eigenvector is stable in time along the same tremor and for all the events. The retrieved mechanism is that of a crack dipping around 20 degrees to the north-northeast and strike of 110 degrees. A plausible interpretative model considers this crack as a fracture in the plug that opens and closes periodically.

This crack acts as a valve through which pulses of gas escape and it is controlled by pressure oscillations in the underlying resonating conduit (Lesage et al., 2006). However it is important to notice that the solutions obtained reveal only that the source mechanisms are predominately isotropic, hence based only on moment tensor solutions and without accurate velocity models we cannot invoke more reliable details on the source processes. As shown by Bean et al. (2008) and De Barros et al. (2011), more reliable results can be obtained using an optimum station distribution which includes stations located close to the source. Therefore more evidence and further studies are needed to provide an exhaustive and detailed model of source mechanism acting at Arenal. These further studies should focus on (i) determining accurate velocity models of volcanoes, (ii) analytically (qualitatively and quantitatively) defining different observable parameters for each proposed source model, and (iii) looking for those observables in extensive datasets of volcanic seismicity in order to validate or falsify those models.

Furthermore, it is important to note that because of the lack of a fully detailed velocity model and an optimum station distribution, the spectral peaks of higher frequency, observed in the tremor spectra, typical of tremor at Arenal, cannot be taken into account in the inversion procedure. For this reason, the retrieved mechanisms are related to the lowest frequency range of the events and present a lowpass filtered representation of the tremor source process.

Conclusion

In this paper we present the results of a moment tensor inversion performed in order to retrieve the source mechanisms generating tremor events at Arenal volcano, Costa Rica. 15 tremor signals with very clear, but emergent, first onsets have been selected from a rich and complex data set, recorded in February 2005. The moment tensor procedure was tested by performing synthetic tests using the whole length and the partial lengths of the signals. Stable and correct results have been obtained when the whole length of the tremor and the correct velocity model have been used, while, when only part of the signal was analysed, the mechanism is correctly reconstructed only for the initial part of the signal.. The moment tensor inversion of the real events was performed in a source location identified through the evaluation of the misfit value using moment tensor grid search. The source is located in a shallow position under the crater summit. Solutions of moment tensor components, with and without single forces, were analysed. The results show a strong stability in time and for all the tremor events with a high isotropic percentage (around 50%) in the solutions.

The source geometry for the tremor signals is represented by a sub-horizontal crack dipping ~20° towards the north-northeast. However more evidence and further studies are needed to obtain better constraints on the source mechanism acting on Arenal volcano. Better station distribution, with some located close to the source and a more reliable velocity model can be the key to obtain such constraints. In the bottom panel, positive gliding is shown. This phenomenon is quite common in Arenal and it is thought to be due to changes in pressure of the fluid in the conduit (Hagerty et al., 2000). As expected due to the differences between velocity models used, the amplitude of the source time function is not correctly retrieved. Panel B: Results of the decomposition for maximum amplitude method for MT (top) and MT+SF (bottom). The blue stars represent the percentage of the isotropic components (85%), while the red triangle and the green crosses represent the percentage of the CLVD (10~%) and DC (10~%) components respectively.

Figure Captions

Introduction

Volcanoes are generally characterised by a wide range of seismic events. This variety is normally due to complexity within these types of environments where different phenomena such as fluid movement, rock fractures, heterogeneity, topography, and near surface velocity structures can generate and affect the seismic response. Volcano-tectonic (VT), long-period (LP) and tremor are some examples of seismic records typical of volcanoes (McNutt, 2005). While the source mechanisms generating VT events are generally well-known, LP and tremor source mechanisms are not fully understood. The generation of these types of events is currently thought to be associated with fluid movements (Chouet, 2003), but a full understanding of these source processes is still very difficult to obtain. Since these types of events may be used as precursory phenomena for volcanic eruptions, a full understanding of their generative mechanisms is of crucial importance to constrain the dynamics of the plumbing system of restless volcanoes. One of the most common tools to derive the seismic source mechanisms is a moment tensor inversion. This approach assumes that the moment tensor components can be represented as a system of equivalent forces that produces the same wavefield as the actual physical processes acting at the source. Moment tensor waveform inversions have been applied to a wide range of volcanoes around the world (Chouet et al., 2003;Nakano et al., 2003;Kumagai et al., 2005;Lokmer et al., 2007;De Barros et al., 20102011). However, while inversions for very long period events (VLP) have been successfully performed (Ohminato et al., 1998;Chouet et al., 2003), the shorter wavelengths of LP events are generally distorted by strong topography and structural heterogeneities. These distortion effects can influence the inversion results, leading to erroneous solutions (Bean et al., 2008). Although the problems related to the inversion of LP and VLP events have been widely studied (Nakano et. al., 2007;Cesca and Dahm, 2007;Bean et al., 2008;Davi et al., 2010), currently there have been no attempts to invert tremor. This is probably due to the difficulty in isolating a clear first onset of the tremor waveforms, the importance of which derives from the idea that the initial part of the signal can be considered to represent a source effect only, while over time the signal becomes more influenced by path effects. Furthermore, the long duration and the complexity of tremor make the inversion procedure, and consequently retrieving the source mechanism, a challenging task. In this article, we present the results of the inversion of 15 tremor events selected from a dataset recorded on Arenal volcano (Costa Rica) during a seismic experiment carried out in 2005. At Arenal, it is possible to isolate the first onset of tremor, hence these events provide a good opportunity to broaden our knowledge about tremor source mechanisms. The source position is determined by performing moment tensor grid searches within a volume of possible source points located under the crater summit. A simultaneous grid search for the source mechanisms and location is performed for each event. The source location is identified through the evaluation of the misfit values where the lowest value is considered the best estimation. Subsequently we defined a common source position evaluating the joint probability density function of the misfit values for all the events. For this common source point, solutions for moment tensor components (MT) and moment tensor components plus single forces (MT+SF) were evaluated. The rationale in choosing As shown later in the text, it is sensible to assume that all the tremor signals were generated at a common source location is based on the fact that it is sensible to assume that all tremor events were generated within the same source, in addition as shown by previous studies on moment tensor inversion, there is a strong trade-off between the source mechanism and source position. The solutions have been analysed using a sliding time window over the retrieved source time function in order to closely analyse the recovered mechanisms. This procedure is tested assessed using synthetic tests in order to evaluate the reliability of the inversion procedure when the sliding time window is used across 5 minute long synthetic signals. The same procedure is then applied to the real data allowing us to investigate possible temporal variations or stability of the source mechanism generating the tremor events in Arenal. Results show temporal stability of the source process generating the selected events. Furthermore, decomposition analysis is performed, which involved analysing the maximum amplitude of the signal considering the polarity of each moment tensor component; these results indicate the presence of a strong and constant isotropic component in the solutions.

Arenal and tremor models

Arenal is a 15 km 3 andesitic strato-volcano located in north western Costa Rica and rises to 1750 m a.s.l., 1100 m elevation with respect to its base (Alvarado and Soto, 2002, figure 1). Since 1968 Arenal is characterized by permanent activity accompanied by gas emissions, Strombolian-Vulcanian explosions, lava and pyroclastic flows. Arenal"s seismicity is mainly characterized by long period events, explosions and tremor. All these kinds of event, including tremor, can be accompanied, or not, by ash emission and sometimes ash can be emitted aseismically (Valade et al., 2011). Tremor events are the most common type of event in Arenal with durations of up to several hours per day. Depending on the frequency content, two different types of tremor, spasmodic and harmonic, can be identified in Arenal, where spasmodic tremor is usually characterized by a wide range of frequencies from 1 to 6 Hz while harmonic tremor spectra contain evenly spaced peaks with fundamental frequency in the range [0.9-2] Hz (figure 2). They both display large temporal variations with gliding phenomena (positive and negative) and jumps in frequency that correspond to sudden decreases following positive gliding (Lesage et al., 2006; figure 3). Different models have been proposed to explain the origin of tremors in Arenal. These models can be divided into three main categories.

According to Benoit and McNutt (1997), tremor events are produced by a fluid-filled resonator (like an organ pipe), where the source fluid is a gas-charged magma. The organ pipe model could be used to explain the evenly spaced peaks of the harmonic tremor. Small variations in the gas content produce strong velocity changes that explain the observed frequency fluctuations. The excitation of resonance is not included in the model.

A second category includes models in which tremors are thought to be generated by the non-linear excitation of a visco-elastic channel caused by different pressure states of a fluid flowing inside the conduit (Julian, 1994;2000). The non-linear excitation is used to model changes in fluid-flow that could account for some of the seismic signals observed at Arenal. These models are based on the assumption that the high velocity fluid is flowing through conduit constrictions. However according to Lesage et al., (2006 and references therein) using a value of a few meters for the conduit diameter and considering the low effusive rate of about [0.3-0.6] m 3 s -1 the resulting magma flow velocity is of the order of a few millimeters per second. The non-linear excitation is used to model changes in fluid-flow that could account for some of the seismic signals observed at Arenal. Therefore the physical conditions that could produce a non-linear excitation are not fulfilled in this volcano (Lesage et al., 2006).

One of the most recent models (Lesage et al., 2006) proposes that tremor events can be generated in a conduit filled with bubbly magma, where the upper part is blocked by solid or very viscous lava acting as a barrier for the gases. The plug is fractured in a way that allows gases to escape. The tremor is generated by repetitive pressure pulses that interact with standing waves in the conduit, a similar process as a clarinet. A detailed analysis of several mechanisms that can generate seismic waves by a fluid flow in a narrow channel suggests that the clarinet model is the most plausible one in volcanoes (Rust et al., 2008). The behaviour of the system is closely related to the properties of the lava plug. These properties are continuously changing with the intrusion of fresh magma into the lava pool and on the way the lava cools down and solidifies to create the plug. The later model is also consistent with joint observations of explosions and tremors by seismometers and Doppler radar (Valade et al., 2011). Their analysis shows large variability of the phenomenon and non-systematic relationship between seismic and radar signals. Tephra emissions were observed during part of the explosions and also during some episodes of tremor and seismic quiescence (Valade et al., 2011).

In this paper, we intend to investigate further and provide additional information on the seismic source mechanisms generating harmonic tremors, in order to broaden our knowledge about the plumbing system.

Data and experimental set-up

In February 2005, a seismic experiment was carried out on Arenal volcano by the University of Costa Rica and the Université de Savoie in France. During the experiment ten broadband seismic stations were deployed around the volcano. Three-components Guralp CMG40T 60-secondsseismometers with mini-Titan recorders were operating in continuous mode at 125 samples per second. Arenal is a dangerous environment due to the strong explosive activity, therefore it is not possible to access the top of the volcano since 2000. The hazardous environment, together with the steep topography and the rich vegetation on the eastern flank of the volcano, strongly influenced the station distribution adopted in the experiment (figure 1). The highest station, ARE1 (at ~750 m a.s.l) is located on the western flank at roughly 2 km from the summit. A pressure sensor and a microphone were also positioned in the same location. The experiment lasted 10 days (from the 11 th to the 21 st of February), during which a wide range of seismic events was recorded (e.g. explosions, LP events, and tremor). The main energy of the different events ranges between 0.8 and 5 Hz. The most common type is tremor, which can vary in duration from up to several hours to days. Often the tremors appear to be triggered by an explosion quake, however the triggering mechanism is not always clear.

We selected 15 harmonic tremor signals with fundamental frequency of 0.8 to 1.5 Hz. The complex behaviour of the volcano makes the analysis of the dataset a difficult task. In addition, due to instrumentation problems, the number of recording stations is not constant for the whole length of the experiment, changing from a minimum of seven to a maximum of ten. For each of the selected tremor events, deconvolution of the instrument response is performed. Filtering, spectrum and spectrogram analysis have been carried out in order to gain a better understanding of the tremor events. Their spectrograms are characterized by spectral lines at different frequency intervals. The evenly spaced overtones are interpreted as resulting of a Dirac comb effect due to a stable repetitive source process (Lesage et al., 2006). Moreover, the spectra are made more complex by site and propagation effects due to heterogeneities in the medium and strong topography (Mora et al., 2001;2006). For example, while the frequency peaks are the same at different stations, strong amplitude variations of the spectral peaks can be observed between close stations. For each of the 15 tremor events, the whole length has been selected. The duration of the tremor is not constant for all the events and it varies from around 100 seconds to more than 500 seconds. Mean and trend have been removed from the dataset. A cosine taper and a bandpass filter (0.8-1.5 Hz) were applied prior to the inversion in order to isolate the fundamental frequency.

Location of tremor events

Generally, the moment tensor inversion procedure is sensitive to a wide range of effects commonly seen in volcanic environments such as those associated with topography, near surface structures and heterogeneities (Davi et al., 2010). In order to reduce the uncertainties intrinsic in the inversion procedure, the number of free parameters must be decreased as much as possible (Bean et al, 2008).

Retrieving the source location could be a way of reducing the free parameters involved in the inversion. The location of the seismic source is crucial to gain a better understanding of the volcanic system and to improve the eruption forecasting (Métaxian et al., 2002). Tremor events are generally characterized by emergent onsets with waveforms affected by scattering effects due to the propagation in through a highly heterogeneous medium and the interaction with the topography. All these characteristics make the use of the classical location methods (such as particle motions or similarities between waveforms) ineffective and imprecise. For this reason, in this paper, the location of the source is estimated through the probabilistic approach (described further in the text) using grid searches performed for the whole data length of 15 tremor events. We created a volume of 4735 possible source points (with a grid step of 15 m), located under the crater summit (480 x 380 x 840 m 3 ), the dimensions of which are have been retrieved by in previous studies of source location (Benoit and McNutt, 1997;Hagerty et al., 2000;Métaxian et al., 2002;2009). For each tremor event, a moment tensor grid search is performed in order to retrieve the waveform misfit values (R) for each point of the source volume and for the solutions with and without single forces.

For the linear inverse problem d=Gm, the misfit between observed (d) and calculated data (Gm) in the least square sense can be expressed by:

    d d Gm d Gm d R T T    (1)
where d , G , m are the data, the Green"s functions and their derivatives, and the moment tensor and single force components, respectively. The superscript "T " indicates the transposed matrix.

The lowest misfit is considered the best estimation.

The single forces can be either considered (the solutions are moment tensor and single forces, MT+SF) or omitted (solution are only moment tensor components, MT). The Green"s functions and their derivative were calculated over the above mentioned source zone of 4735 points using 3D full wavefield numerical simulations (O"Brien and Bean, 2004). The topography and 1D velocity model are included in the simulations. The one-dimensional velocity was retrieved from sounding using the spatial autocorrelation method (SPAC) and seismic refraction experiments carried out on Arenal in 1997 (Mora et al., 2006). The numerical domain consists of a 13 x 11 x 6 km 3 where the topography is derived from the Digital Elevation Model (DEM) of the volcano using a spatial grid step of 15 m. The Green"s functions and their derivative are calculated using the Reciprocity theorem (e.g. Aki and Richards, 2002). The method and the inversion techniques are described in Davi et al. (2010). The result of this calculation is a set of 15 second long Green's functions with a frequency of up 5 Hz.

We analysed the time delays between the signals recorded at different stations and found that these delays are equal for each event. These results suggest a common source location for all the events.

With this assumption we performed grid searches for each tremor event, results of which indicate a point of minimum misfit for both solutions (moment tensor, MT and moment tensor plus single forces, MT+SF, respectively). To find a common source point for all the tremor events, we evaluate the joint probability density function P(dj) i.e. the likelihood that the common source is located at a certain position of the considered volume (Menke, 1984):

           N = i j i j ) (d R d P 1 2 1 exp (2)
where d j denotes the j-th source point, and N is the number of events. In order to find a common source position, inversions for MT and MT+SF for each of the 15 tremor events have been calculated, and eq. ( 2) is then applied to find the joint probability density function P(d j ). By evaluating the maximum of the joint probability density function, we were able to locate the common source position in an area of the source volume located under the crater summit for MT solution (the minimum misfit value for MT+SF solution is located nearby). Due to the similarity of the MT and MT+SF solution, and since Davi et al. (2010) stated that in the case of Arenal and with this station configuration the MT solution is more stable, we decided to fix the source location in the position retrieved by the grid search for the MT solution (figure 4, panel A). Figure 4, panel B

shows the area of the volume with a 90% percentage of probability of finding the real location of the source. Since this area is large (it has an extension of 100 m in each direction) we can confirm the reliability of our source location. The source is located 100 meters below the crater summit. The results of the grid search indicate that the retrieved position is in agreement with the results of Métaxian et al. (2002;2009). According to these authors, the source position of the tremor is located in a small area centred under the crater summit.

Synthetic tests

The inversion of tremor can be a sensitive procedure; in fact the complexity and the long duration of these types of events, along with errors in the velocity model, can easily influence the inversion results. For these reasons particular care must be taken in the interpretation of the results. Davi et al.

(2010) performed some synthetic tests in order to verify the consistency and limitations related to the inversion procedure on Arenal volcano. The tests were performed using the best estimate of the velocity model and a station configuration which mapped the real location of ten stations deployed around Arenal volcano in 2005. The results highlighted how spurious single forces are easily generated under conditions common on volcanoes, such as noisy data, imperfect knowledge of the velocity model, and mislocated source positions. The same results are considered valid for the inversion of tremor events. These phenomena, common of volcanic environments, can influence the correct retrieval of the source time function and lead to erroneous mechanisms. This information was used in the interpretation of the results of an inversion for an explosive event recorded on Arenal in 2005. Since tremor events can be seen as the product of a series of repetitive source time functions and, thus, differing only in their longer duration, the same results are considered valid for the inversion of the tremor signals. Hence, only the moment tensor components for moment tensor (MT) and moment tensor plus single forces (MT+SF) solutions were analysed while the single forces were neglected (see also De Barros et al., 2011).

Compared to the previous study (Davi et al, 2010), an additional difficulty in inverting tremor is the length of the signal. It is computationally very expensive to invert signals that last for hourswith duration of hours. Hence by using synthetic data we can derive the best strategy to find the source mechanism of long lasting signals. For these reasons, we tested the inversion procedure designing two different types of synthetic tests in which (i) we inverted for the whole length of a synthetic tremor in order to test how the inversion procedure performs when long windows of data are utilized, and (ii) performing the inversion using only the initial part of the signal. The solutions were analyzed using a short sliding window (of 32 s) in order to verify the stability or variability over time of the mechanism generating the tremor events. This analysis can be made using the maximum amplitude of the signal (denoted MAS in this article), taking into account the polarity of each component (Chouet et al., 2005). An additional test is performed to investigate the sensitivity of the inversion to a poorly detailed velocity model.

From the convolution between a Dirac Comb comb and a Ricker wavelet function with a central frequency of 2Hz, we created a 5 minute long synthetic tremor signal with the same spectral characteristics as the recorded signals. The source position is in the location point found through the grid searches. We simulated two source mechanisms (i) a pure volumetric (M xx =M yy =M zz =1;

M xy =M xz =M yz =0) and (ii) a compensate linear vector dipole (CLVD) mechanism (M xx =M yy =1;

M zz =-2, M xy =M xz =M yz =0).

Inversion of the full length of the signal

In the first test an explosive mechanism was simulated using a 5 minute signal. The longer the duration of the tremor, the larger is the required computational time. However the solutions for the entire length appear to be very stable with a misfit value approaching zero (figure 5). With a perfect velocity modelT the source time function and mechanism are perfectly retrieved. We selected a 32 second long time window and we slide this window through the solution. From these results we can conclude that, despite using Green's functions which are simpler than the true model, the relative percentage of ISO seems to be correctly retrieved. On the contrary the deviatoric part of the solution shows more inaccuracies. However the relative percentages of ISO, CLVD and DC remain constant along the solution which further validate the use of our decomposition method in order to assess the stability or variability of the source mechanism.

Inversion of partial length of the signal

The situation changes when only the initial part of the signal is used in the inversion. As the whole length of the signal cannot be always inverted because of the computation cost, the initial part of the signal is chosen because this is considered to be less affected by scattering and uncertainties due to for CASE IV).. In addition, spurious single forces and off-diagonal components are generated in the interval in which the retrieved solution is no longer correct (figure 8 shows results for CASE IV).

This suggests that in the inversion procedure, the errors are easily mapped in the single forces and the off-diagonal components.

The observed inability to obtain a correct solution for the final part of the incomplete signal can be easily explained if the tremor signal is seen as a set of consecutive sources (like a Dirac comb).

Since each of these consecutive sources has a propagation effect of the same length as our Green"s functions, only those sources for which the whole length of the propagation effect is included in the signal can be successfully retrieved. In other words, the retrieved source-time function will be correct within a window which starts at the beginning of the signal and whose length is equal to the difference between the length of the signal and the length of Green"s functions.

In summary, in cases when the exact velocity model and the entire length of the signal are is used, the correct source mechanism and source time function are retrieved for the whole length of the signal. This allows us to analyse such a solution within multiple sliding time windows in order to verify the stability or the variability of the mechanism generating the tremor. Hence the use of the whole length of the signal will allows us to detect, if present, variations in the source mechanisms and to analyse the influence of propagation effects through the signals. On the contrary, when inverting for a part of the signal, only the beginning part (length of the data minus length of the Green's functions) of the retrieved source-time function can be considered as a correct solution.

Therefore this approach can be useful for inverting very long tremors, when the full length of the signal cannot be included in the inversion procedure.

Tremor inversion results and discussion

We selected 15 real tremor events recorded in Arenal volcano and performed a moment tensor inversion to retrieve the source mechanisms generating these events. The waveform inversion of the tremor is performed in the frequency domain for moment tensor solutions with and without single forces. For each of the selected events, the whole length of the tremor has been used to perform the inversion; the shortest events have durations of 100 s while the longest ones last for 600 s. All the selected tremor signals show a clear but emergent onset that can be easily isolated. From the analysis of the data, it is clearly noticeable that the waveforms differ from one event to another in length and amplitude. However the solutions of the inversions for each tremor show very similar

patterns. An example of the moment tensor solution is shown in figure 9 in which the obtained MT+SF and MT solutions are plotted for the event recorded on February 13 th at 13:00. The fit between calculated and observed data is illustrated in figure 10. The fit (R= 0.6) shows a reasonable match for nearly all the stations which is indicative of the stability and accuracy of our inversion solution. In some cases the difference in amplitude between the observed and calculated data can be originated by local site effects (see station 4). These results are similar for the 15 events. However, in the MT+SF solution large single forces, especially in the z direction, are generated. In addition, large amplitude off-diagonal components in both MT and MT+SF solutions are retrieved. From the results obtained by the synthetic tests, (Davi et al., 2010), we highlighted how, in the case of Arenal, the errors caused by the poorly constrained velocity model, and the non-ideal station distribution, easily generated spurious off-diagonal components and single forces (Davi et al., 2010). Hence, we cannot easily distinguish, when inverting the real dataset, if, and to what extent, these components are real or spurious. For these reasons we excluded the single forces from the decomposition procedure. The same conclusions were reached by De Barros et al. (2011) for Mt. Etna volcano, Italy. We selected a constant window of 32 seconds and performed the decomposition sliding the window over the whole length of the solution for each event. For each window the percentage of isotropic, compensate linear vector dipole and double couple are calculated in order to retrieve the source mechanism generating the events. Results of the decomposition are shown in figure 11 for MT and figure 12 for MT+SF solutions (shown here for comparison with the MT solution). For each tremor event a stable isotropic percentage (ISO%) is found. The ISO% varies from 40% to 60% with some sporadic examples at 30%. Therefore the ISO% shows a small variation with a standard deviation in the order of the 5% from the mean value. These results are common for the moment tensor components of both solutions, with and without including the single forces. The CLVD and DC components show a more scattered pattern, especially for the solution including the single forces. The CLVD component shows a standard deviation in the order of 12% while it reaches 17% in the case of the double couple components. The same results are obtained (not here) when shorter and longer (e.g. 16 seconds and 40 seconds, respectively) decomposing time are utilized. Results were also tested with an overlapping window (length of 24 s overlapping every 12 s). No significant variations were found. In addition theThe stability of the retrieved source mechanisms indicatehighlights thate reliability of our solutions which do not seem to be are not strongly affected by path effects.

As shown by the synthetic tests, spurious off-diagonal components are easily generatedincorrect source mechanisms can be retrieved . In addition, complex decomposition methods do not always show a stable behaviour when there are significant uncertainties in the data. Since the level of uncertainty increases (due to poorly detailed velocity models, heterogeneities, network configuration, or source location errors) in the case of real waveform inversions, the percentage of these spurious off-diagonal componentsincorrect source mechanisms can increase. Since the decomposition procedure of the deviatoric part of moment tensor (i) is not uniquely defined and (ii) can be very sensitive to the noise contaminating MT solutions, the interpretation of the mechanism can be ambiguousnly the volumetric part of a source time function is reliably retrieved while further studies are needed for a complete understanding of the deviatoric part (CLVD%+DC%) of moment tensor solutions. This is confirmed by the results of the decomposition of our real data. Furthermore, the outcome of the decomposition confirms that, in the case of Arenal volcano, with this station configuration and with the velocity model we employed, the solutions without the single forces are more stable than the MT+SF solutions. The same conclusion was reached previously (Davi et al, 2010) when the inversion procedure was tested for explosive events. For this reason we will base the interpretation of our results on the evidence shown by the MT solutions.

The decomposition method indicates solutions with 50 ± 5% of isotropic component, 35 ± 12% of CLVD and 1715 ± 17% of DC for all the events. The overall trend shows a solution which is very stable in time and from one event to another. This confirms that the same non-destructive source processes generate the tremor and that the mechanisms do not vary in time. The diagonal components are, for each event, in a ratio of roughly [1 1 2]. Since the deviatoric part is not perfectly solved and our solutions show large uncertainties in the DC percentage of the 17%, these results could be interpreted as the mechanism of a crack that is opening in a direction nonperpendicular to the crack plane. The analysis of the eigenvectors (figure 13) shows that a subhorizontal crack with constant orientation is found for the 15 tremor events. The orientation of the eigenvector is stable in time along the same tremor and for all the events. The retrieved mechanism is that of a crack dipping around 20 degrees to the north-northeast and strike of 110 degrees. A plausible interpretative model considers this crack as a fracture in the plug that opens and closes periodically. This crack acts as a valve through which pulses of gas escape and it is controlled by pressure oscillations in the underlying resonating conduit (Lesage et al., 2006) Furthermore, it is important to note that because of the lack of a fully detailed velocity model and an optimum station distribution, the spectral peaks of higher frequency, observed in the tremor spectra, typical of tremor at Arenal, cannot be taken into account in the inversion procedure. For this reason, the retrieved mechanisms are related to the lowest frequency range of the events and present a lowpass filtered representation of the tremor source process.

Conclusion

In this paper we present the results of a moment tensor inversion performed in order to retrieve the source mechanisms generating tremor events at Arenal volcano, Costa Rica. 15 tremor signals with very clear, but emergent, first onsets have been selected from a rich and complex data set, recorded in February 2005. The moment tensor procedure was tested by performing synthetic tests using the whole length and the partial lengths of the signals. Stable and correct results have been obtained

when the whole length of the tremor and the correct velocity model have been used, while, when only part of the signal was analysed, the mechanism is correctly reconstructed only for the initial part of the signal.spurious single forces and off-diagonal componentsincorrect source mechanisms were retrieved in the final part of the source time function, when only part of the signal was analyzed. The moment tensor inversion of the real events was performed in a source location identified through the evaluation of the misfit value using moment tensor grid search. The source is located in a shallow position under the crater summit. Solutions of moment tensor components, with and without single forces, were analysed. The results show a strong stability in time and for all the tremor events with a high isotropic percentage (around 50%) in the solutions. The source geometry for the tremor signals is represented by a sub-horizontal crack dipping ~20° towards the north-northeast. However more evidence and further studies are needed to obtain better constraints on the source mechanism acting on Arenal volcano. Better station distribution, with some located close to the source and a more reliable velocity model can be the key to obtain such constraints. In the bottom panel, positive gliding is shown. This phenomenon is quite common in Arenal and it is thought to be due to changes in pressure of the fluid in the conduit (Hagerty et al., 2000). The distances are expressed in meters.

  uncertainties due to the unknown velocity model. The inversion has been performed with different intervals of the signal; for CASE I we generated 20 s of synthetic signal, in CASE II, 25 s, in CASE III, 30 s and in CASE IV (shown in figure7), 35 s. A simpler source time function is used to show the results obtained by the inversion. For each solution (MT and MT+SF), the initial part is well retrieved in terms of both source time function and mechanism, however at a certain point the solutions becomes unstable, showing spurious single forces and incorrect source mechanism. The same results are obtained in the four cases, with the time window in which the source time function and mechanism are correctly retrieved increasing with the length of the signal (the first 10 s for CASE I, 15 s for CASE II, 20 s for CASE III and 25 s for CASE IV, were correctly retrieved, respectively). The results from the decomposition show that for the initial part the explosive ratio of [1 1 1], a 98% of ISO and around 2% of CLVD and DC is retrieved. The final part of the retrieved solution, however, does not show the correct ratio of [1 1 1] (figure8shows results for CASE IV).

Figure 1 .

 1 Figure 1. Network configuration deployed in 2005 during a seismic experiment carried out at

Figure 2 .

 2 Figure 2. Examples of a harmonic tremor recorded in 2005 on February 14 th at 8:00. The panels

Figure 3 .

 3 Figure 3. Example of a spasmodic tremor recorded in Arenal in 2005 on February 14 th at 13:00. The

Figure 4 .

 4 Figure 4. Panel A: Normalized joint probability density function of the source position for the 15

Figure 5 .

 5 Figure 5. Moment tensor inversion solution (amplitude expressed in 10 -5 Nm) for MT+SF (left

Figure 6 .

 6 Figure 6. Decomposition results of the solution without single forces. The MT+SF solution, not

Figure 7 .

 7 Figure 7. Panel A: Moment tensor inversion solution for MT+SF (left) and MT (right) when

Figure 8 .

 8 Figure 8. Moment tensor inversion solution (amplitude expressed in 10 -5 Nm) for MT+SF (left

Figure 9 .

 9 Figure 9. Decomposition results of the solution with (top panel) and without the single forces

Figure 10 .

 10 Figure 10. Example of a moment tensor solution of a tremor event recorded on February 13 th at

Figure 11 .

 11 Figure 11. Reconstructed (red line) and observed seismograms (blue line) are compared for the waveform inversion of the tremor event that occurred in 2005 on February 13 th at 13:00 for solution

Figure 12 .

 12 Figure 12. Results of the decomposition for the tremor events inversion (MT solution) are shown in

Figure 13 .

 13 Figure 13. Same as figure 11 for the MT+SF solution. Misfit values are indicated in red on the left

  For each time window we decompose the solutions of MT and MT+SF components. The decomposition results of the pure volumetric mechanism show a [1 1 1] ratio of the moment tensor diagonal components (M xx , M yy , M zz ) in every interval for both solutions, MT and MT+SF. The percentage of isotropic components (% ISO) is roughly 100%, while the double couple (DC) and compensate linear vector dipole (CLVD) percentages are very close to zero. Results are shown in figure 6. The same test (not shown here) was performed using the pure CLVD mechanism [1 1 -2] and again the results of the decomposition indicate a 100% of CLVD and a 0% of ISO and DC components. When perfect conditions (correct velocity model and the whole length of signal) are used, the mechanism and the source time function are perfectly retrieved.A second test was performed in order to assess the role played by poorly detailed velocity model when long and complex waveforms are used in the inversion procedure. Due to lateral heterogeneities and potential scatters the real velocity model is commonly more complex than what we use for modelling. Consequently, to resemble a real-like situation in order to test the influence of the velocity model on the inversion results we used a set of Green's functions calculated for a simpler velocity model (e.g. a homogeneous velocity model), to invert a signal calculated for a more realistic model (the layered model derived byMora et al., 2006). We simulated a geometry of an horizontal crack and we performed the inversion using the 5 minutes synthetic signals. Results of the inversion for MT and MT+SF are shown in figure7 (panel A). The solutions have been subsequently decomposed using again a window of 32 seconds. Figure 7 (panel B) shows the decomposition results. The solution without including the single forces (MT) appears to be more Formatted: Font color: Auto reliable showing a dominant isotropic component within ~65%, which remains constant along the solution.

  the unknown velocity model. The inversion has been performed with different intervals of the signal; for CASE I we generated 20 s of synthetic signal, in CASE II, 27 25 s, in CASE III, 30 s and in CASE IV, 35 s (shown in figure 7), 35 s. A simpler source time function is used to show the results obtained by the inversion. For each solution (MT and MT+SF), the initial part is well retrieved in terms of both source time function and mechanism, however at a certain point the solutions becomes unstable., showing sSpurious single forces and off-diagonal components are generatedincorrect source mechanism. The same results are obtained for in the four cases, with the difference being when a longer signal is used in the inversion, the time window in which the source time function and mechanism are correctly retrieved, increasing with the length of the signales (the first 10 s for CASE I, 15 s for CASE II, 20 s for CASE III and 25 s for CASE IV, were correctly retrieved, respectively). The results from the decomposition show that for the initial part the explosive ratio of [1 1 1], a 98% of ISO and around 2% of CLVD and DC is retrieved. The final part of the retrieved solution, however, does not show the correct ratio of [1 1 1] (figure 8 shows results

Figure 1 .

 1 Figure 1. Network configuration deployed in 2005 during a seismic experiment carried out at Arenal. The black diamonds represent the locations of the 60-second Guralp CMG40T seismometers. The red star indicates the source location of the tremor signals analysed in this study. Maximum elevation shown is 1600 m (100 m between contour lines). In the bottom left panel Arenal volcano and its location in north-western Costa Rica are shown.

Figure 2 .

 2 Figure 2. Examples of a harmonic tremor recorded in 2005 on February 14 th at 8:00. The panels show, from the top to bottom respectively, the waveform, the spectrum and the spectrogram. The

Figure 3 .

 3 Figure 3. Example of a spasmodic tremor recorded in Arenal in 2005 on February 14 th at 13:00. The panels show, from the top to bottom respectively, the waveform, the spectrum and the spectrogram.

Figure 4 .

 4 Figure 4. Panel A: Normalized joint probability density function of the source position for the 15 tremor events for moment tensor components without single forces (MT). The dark star indicates the chosen source position. Panel B shows the area where the real source is more likely located (probability of 90% with maximum distance between the sides of the surface in the order of 100 m).

Figure 5 .

 5 Figure 5. Moment tensor inversion solution (amplitude expressed in 10 -5 Nm) for MT+SF (left panel) and MT (right panel) when 5 minute long synthetic signal is used in the inversion. This signal is obtained by the convolution between a Dirac Combcomb, and a Ricker wavelet function. The original source time function (shown in the right top panel) and the explosive mechanism are perfectly retrieved by the inversion.

Figure 6 .

 6 Figure 6. Decomposition results of the solution without single forces. The MT+SF solution, not shown here, is identical. Results of the decomposition for maximum amplitude method. The blue stars represent the percentage of the isotropic components (100%), while the red triangle and the green crosses represent the percentage of the CLVD (~0%) and DC (0%) components respectively.

Figure 7 .

 7 Figure 7. Panel A: Moment tensor inversion solution for MT+SF (left) and MT (right) when complex synthetic signals and homogeneous Green"s functions are used in the inversion procedure.

Figure 78 .

 78 Figure 78. Moment tensor inversion solution (amplitude expressed in 10 -5 Nm) for MT+SF (left panel) and MT (right panel) when 35 seconds of the signal are included in the inversion. The isotropic mechanism is perfectly retrieved only for the beginning part of the solution (~20 seconds over a total length of 35 seconds), while in the final part spurious single forces and off-diagonal components are generated.

Figure 89 .

 89 Figure 89. Decomposition results of the solution with (top panel) and without the single forces (bottom panel) included in the inversion. The red triangles, the blue stars and the green crosses represent the percentage of the CLVD, of ISO and DC components respectively. The right percentage, e.g. 100% ISO and 0% of CLVD and DC (pure isotropic mechanism) are obtained only for the initial part while the relative ratio changes after the first 20 second of signal. This means that the mechanism is no longer correctly retrieved by the inversion.

Figure 910 .

 910 Figure 910. Example of a moment tensor solution of a tremor event recorded on February 13 th at 13:00 Panel A: Moment tensor inversion solution for MT+SF (left) and MT (right) are shown. (see in Panel B: where from the top to bottom the waveform, the spectrum and the spectrogram are shown, respectively). Panel A: Moment tensor inversion solution for MT+SF (left panel) and MT (right panel) are shown.

Figure 1011 .

 1011 Figure 1011. Reconstructed (red line) and observed seismograms (blue line) are compared for the

Figure 1112 .

 1112 Figure 1112. Results of the decomposition for the tremor events inversion (MT solution) are shown in panel A. Misfit values (60%) are indicated in red on the top right-hand side. An overall constant percentage of the ISO components can be seen for each event. The top, medium and bottom plots in panel B show the percentage of ISO (blue stars), CLVD (red triangles) and DC (green crosses), respectively. Horizontal lines indicate mean values and mean values +/-one standard deviation for all the time intervals of all the tremor. The CLVD and DC components are more scattered and more variable than the ISO percentage. The x axis indicates the number of decomposing window for each event, increasing progressively in time.

Figure 1213 .

 1213 Figure 1213. Same as figure 11 for the MT+SF solution. Misfit values are indicated in red on the left of each plot of Panel A. The figure is plotted to give a comparison with the MT solution

  . However it is important to notice that the solutions obtained reveal only that the source mechanisms are predominately isotropic, hence based only on moment tensor solutions and without accurate velocity models we cannot invoke more reliable details on the source processes. As shown byBean et al. (2008) and DeBarros et al. (2011), more reliable results can be obtained using an optimum station distribution which includes stations located close to the source. Therefore more evidence and further studies are needed to provide an exhaustive and detailed model of source mechanism acting at Arenal. These further studies should focus on (i) determining accurate velocity models of volcanoes, (ii) analytically (qualitatively and quantitatively) defining different observable parameters for each proposed source model, and (iii) looking for those observables in extensive datasets of volcanic seismicity in order to validate or falsify those models.

Acknowledgements

This work has been funded by Science Foundation Ireland (SFI). The authors wish to acknowledge the Irish Centre for High-End Computing (ICHEC) for providing computational facilities. The fieldwork was partly supported by the European Commission, 6 th Framework Project -"VOLUME", Contract No. 018471, INSU-CNRS (ACI Risques naturels et changements climatiques), Université de Savoie (BQR B2005-09), and projects n° 113-A6-503 and 113-A7-511 from Universidad de Costa Rica and by the volcanological studies of the Instituto Costarricense de Electricidad. We thank the staff of Escuela Centroamericana de Geología, Universidad de Costa Rica, Instituto Costarricense de Electricidad for their efficient logistical support. Acknowledgements This work has been funded by Science Foundation Ireland (SFI). The authors wish to acknowledge the Irish Centre for High-End Computing (ICHEC) for providing computational facilities. The fieldwork was partly supported by the European Commission, 6 th Framework Project -"VOLUME", Contract No. 018471, INSU-CNRS (ACI Risques naturels et changements climatiques), Université de Savoie (BQR B2005-09), and projects n° 113-A6-503 and 113-A7-511 from Universidad de Costa Rica and by the volcanological studies of the Instituto Costarricense de Electricidad. We thank the staff of Escuela Centroamericana de Geología, Universidad de Costa Rica, Instituto Costarricense de Electricidad for their efficient logistical support.