A brittle failure model for long‐period seismic events recorded at Turrialba Volcano, Costa Rica

A temporary seismic network, consisting of 23 broadband and six short‐period stations, was installed in a dense network at Turrialba Volcano, Costa Rica, between 8 March and 4 May 2011. During this time 513 long‐period (LP) events were observed. Due to their pulse‐like waveforms, the hypothesis that the events are generated by a slow‐failure mechanism, based on a recent new model by Bean et al. (2014), is tested. A significant number (107) of the LPs are jointly inverted for their source locations and mechanisms, using full‐waveform moment tensor inversion. The locations are mostly shallow, with depths < 800 m below the active Southwest Crater. The results of the decompositions of the obtained moment tensor solutions show complex source mechanisms, composed of high proportions of isotropic and low, but seemingly significant, proportions of compensated linear vector dipole and double‐couple components. It is demonstrated that this can be explained as mode I tensile fracturing with a strong shear component. The source mechanism is further investigated by exploring scaling laws within the data. The LPs recorded follow relationships very similar to those of conventional earthquakes, exhibiting frequency‐magnitude and corner frequency versus magnitude relationships that can be explained by brittle failure. All of these observations indicate that a slow‐failure source model can successfully describe the generation of short‐duration LP events at Turrialba Volcano.


Introduction
Turrialba Volcano (3340 m above sea level (asl)) is a stratovolcano located at the eastern edge of the Cordillera Volcánica Central, Costa Rica, Central America, 35 km ENE of the capital San José (Figure 1a). Activity at the volcano has increased in recent years, with intensifying fumarolic gas discharges, especially in the active Southwest Crater; summit inflation detected from 2005 to 2007; increases in seismicity; and recent minor phreatic explosions [Campion et al., 2012;Martini et al., 2010;Soto and Mora, 2013;Tassi et al., 2004;Vaselli et al., 2010]. Monitoring the volcano is therefore extremely important, as its close proximity to San José means that a large eruption has the potential to adversely affect 1.5 million people [Soto et al., 2010].
A previous study was carried out into the source processes of long-period (LP) seismic events recorded at Turrialba Volcano in 2009 [Eyre et al., 2013], where LPs were inverted for their mechanism using full-waveform moment tensor inversion. Similar to studies conducted at other volcanoes, a crack source mechanism was obtained, in the case of Turrialba dipping shallowly toward the southwest. For the 2009 data set this was interpreted, drawing on the Chouet [1986] model of LP sources, as resonance within cracks in the hydrothermal system that follow weaknesses in the volcano edifice or as hydrothermal fluids "pulsing" through the cracks. However, the authors state that a deformation-based source process could not be ruled out, as this was postulated as the cause of similar events observed on Mount Etna [De Barros et al., 2011].
volcano. Enhanced coverage of the summit area has been demonstrated in previous studies to improve the moment tensor solutions [Bean et al., 2008]; [De Barros et al., 2011]. As well as improved solutions, this is also a more in-depth study than that of Eyre et al. [2013], with moment tensor inversion applied to over 100 events.
During the experiment, Turrialba Volcano was undergoing large-scale degassing, with up to~3000 t of SO 2 per day expelled into the atmosphere [Campion et al., 2012]. The majority of this gas was released from a large fumarole located high on the southwest inner wall of the Southwest Crater, which opened Journal of Geophysical Research: Solid Earth during a minor phreatic eruption in January 2010. Twenty broadband (five 60 s and three 30 s) and six short-period (1 s) seismometers were installed in a temporary network between 8 March and 4 May 2011, in addition to a three-station permanent broadband network (stations VTCE, VTCG, and VTUN). Figure 1b shows a map of the network. The nonactive craters (the Northeast Crater and part of the Central Crater) allow for a dense deployment of stations in the summit region of the volcano. The network consisted of three seismic arrays and 13 stand-alone stations. Seismic arrays were located at TOMA, PIGA, and BABO. TOMA consisted of the six short-period stations in the shape of a cross, with the long axis pointing to the active Central and Southwest Craters and an interstation spacing of approximately 60 m. The arrays BABO and PIGA consisted of five stations each, in a "T" shape with their long axes pointing to the summit region. All stations had a sampling rate of 100 Hz. The station WCR1 was located on the rim of the active Southwest Crater, within 100 m of the main fumarole, and WCR2 is also located on the crater rim. As the LPs in 2009 were located in this region at shallow depths [Eyre et al., 2013], it is likely both stations were very close to the source of the seismicity. However, the close proximity to the strong gas discharge caused these stations to record strong continuous energy in the 5 Hz to 40 Hz range.
There are currently two main categories of source models available for LP events: (1) the resonating crack/conduit model [Chouet, 1986[Chouet, , 1988Neuberg et al., 2006] and (2) a very recently proposed slow brittle rupture model [Bean et al., 2014]. In the first model, the lack of high-frequency content of the events is explained by relating the source processes causing the events to fluid processes. A fluid-filled cavity is excited into resonance, inducing slow dispersive waves (crack waves) which propagate along the cavity interface, interfering constructively to produce the LPs observed. By definition this process leads to resonating seismograms. The second model attributes the low corner frequencies of the events to slow failure of unconsolidated volcanic material at shallow depths on volcanoes caused by the low internal friction angles of the material. Rupture simulations show that this process leads to pulse-like seismograms when not contaminated by path effects.
Examples of LP events recorded in 2011 and their corresponding spectra can be seen in Figure 2, recorded on station WCR1, low-pass filtered at 5 Hz to remove the 3 Hz in red as used in the moment tensor inversion) and the corresponding frequency spectra. The events are low-pass filtered below 5 Hz to remove the strong noise above this frequency believed to be caused by the strong degassing. strong noise due to degassing above this frequency. The main energy of the events is generally from 0.3 to 1.3 Hz. The events exhibit pulse-like waveforms, and therefore, the brittle failure model appears to be more likely to fit the data, as there is no resonance in the waveform which would be expected according to the cavity resonance model. As the LP waveforms on Turrialba Volcano are nonresonating, by analyzing the waveforms, source locations, source mechanisms, and the scaling laws of these events, this study aims to determine if the slow-failure model can reasonably explain the LP events observed on Turrialba Volcano.

LP Event Identification
There were 513 LP events recorded between 8 March and 4 May 2011. These LPs are identified using the STA/LTA algorithm described in Eyre et al. [2013] for station TCCR, where the STA/LTA method is implemented for both a low frequency band (0.2-4.0 Hz) and a high frequency band (6.0-10.0 Hz), and the events are only picked when the STA/LTA value in the low frequency band is 3 times that of the high frequency band. An average of nine events per day is detected.
LPs are filtered in a frequency band of 0.3-1.3 Hz (the frequency band in which they have most energy) and then correlated in an attempt to group them into families of similar events using a closed cluster technique [Eyre et al., 2013;Saccorotti et al., 2007] for 10 s windows of the Z component at summit station TCCR. No families of more than five events with a correlation coefficient greater than 0.75 could be found, suggesting that these events are not produced by a repetitive source. Hence, the events are treated individually.

Method
Following previous studies by Chouet et al. [2003], Kumagai et al. [2005Kumagai et al. [ , 2010, Nakano and Kumagai [2005], Ohminato et al. [1998], and Eyre et al. [2013], simultaneous full-waveform moment tensor inversions to calculate the source mechanisms and grid searches for source locations are performed. The method is implemented identically to Eyre et al. [2013], with the exception of the velocity model. In the frequency domain, the nth component of the displacement u, recorded at position x and frequency ω, can be written as follows: where M pq is the force couple in the direction pq, F p is the single force acting in the direction p, and G np and G np,q are the nth components of the Green's functions and their spatial derivatives with respect to the source coordinates, respectively, generated by the single force F p and the moment M pq , respectively. Green's functions are the displacement responses recorded at the receivers when an impulse force function is applied at the source position in an elastic Earth (the medium response). The summation convention applies. Equation (1) is linear, and in the frequency domain can be written in matrix form where d is the data matrix, G contains the Green's functions, and m contains the moment tensor and single forces components. As the LP events are recorded in the near field, P and S waves are intertwined. The near-field effect is accurately taken into account in our Green's functions calculations, allowing for a full-waveform inversion. The quality of the inversion results is tested through the evaluation of the misfit between the calculated and the observed data, herein referred to as the residual R: where W is a diagonal matrix of weights for the quality of the data (in this study the data quality is good for all stations and components, and therefore, all of the weights are fixed to 1: the inversion intrinsically preferentially weights stations closer to the source). The inversion is carried out in the 0.3-1.3 Hz frequency band used for the correlations (section 2) for 12 s time windows, utilizing all components at all available stations.
In order to locate the events, the inversion is implemented for a grid of 2536 possible source points in a 4000 m by 4000 m by 2000 m volume below the summit region of the volcano Kumagai et al., 2005Kumagai et al., , 2010 and 200 m farther from the summit. The number of possible source positions is much greater than the number of stations used in the inversion (15 stand-alone broadband stations are available for the inversion), so the Reciprocity Theorem [Aki and Richards, 2002] is used when computing the Green's functions between the sources and receivers. The most probable source location is defined by the grid point with the smallest misfit between observed and reconstructed data.
The Green's functions are calculated using a 3-D elastic lattice method to simulate full-wavefield seismic wave propagation [O'Brien and Bean, 2004]. The model includes topography which has been shown to be important for an accurate inversion in volcanic environments [Ripperger et al., 2003]. The shallow velocity structure has also been shown to be important in inversions [Bean et al., 2008]. As no model is currently available for Turrialba Volcano [Eyre et al., 2013], the shallow velocity structure is investigated using surface wave dispersion analysis using the three available arrays (see section 1). The technique used is the multichannel analysis of surface waves (MASW) method [Wathelet et al., 2004]. However, this analysis is very limited due to the small aperture of the arrays (the maximum station spacing is 180 m due to the steep topography of the volcano). Although the results are confined to the very near surface, analysis shows that a low-velocity layer exists at Turrialba Volcano as has been found on many other volcanoes [Bean et al., 2008, and references therein]. Hence, the numerical model used for simulating the Green's functions includes a low-velocity layer of 1039 m/s (P wave) and 600 m/s (S wave), tracking the topography. However, as the surface wave inversion is not well constrained, and we do not want to include an unrealistic poorly constrained step in the model which would act as an unrealistic waveguide, we choose the smoothest model that fits the results. Therefore, the velocities within the layer are gradated to 100 m depth, beneath which the model is homogeneous with a P wave velocity of 2300 m/s and an S wave velocity of 1328 m/s. Density is related via the equation ρ = 1700 + 0.2v p , where v p is the P wave velocity [Gardner et al., 1974;O'Brien and Bean, 2009].
As the inversion is performed in the frequency domain without any assumption on the source-time history, an inverse FFT of the solution m from the inversion of equation (2) leads to a source-time function for every component of the moment tensor and single forces. In order to obtain a common source-time function, the moment tensor solutions are decomposed into their principal components, using the method of Vasco [1989]. This method involves the singular value decomposition of the obtained six time-dependent moment tensor components. A common source-time function for all six moment tensor components and its contribution to each are obtained, thus giving a source-time history of the source process and its mechanism. The eigenvalues of the "static" moment tensor give the pressure and tension values along the principal axes (eigenvectors) from which the source mechanism can be reconstructed. From this, the solution can be decomposed into the percentage of isotropic, compensated linear vector dipole (CLVD) and double-couple source mechanisms [Vavryčuk, 2001]. Hence, for example, a tensile crack source (often calculated as the most likely mechanism for LPs at volcanoes [e.g., Cusano et al., 2008;De Barros et al., 2011]) produces a λ:λ:(λ + 2 μ) ratio of moment tensor eigenvalues (where λ and μ are the Lamé parameters) [Aki and Richards, 2002], which is 1:1:3 for a Poisson's ratio of 0.25 and 1:1:2 for a Poisson's ratio of 0.33.
Inversions are carried out both excluding and including single forces in the solutions. Single forces are physically realistic in volcanic environments (for example, due to mass transfer, drag forces, or volcanic jets) [Chouet, 2003;Takei and Kumazawa, 1994]. However, De Barros et al. [2011 state that errors in modeling the velocity structure mainly contaminate the single forces, leading to spurious single forces, while the moment tensor solution is still correct. In particular, De  showed using radiation patterns that strong spurious forces, of the order of those usually derived from the moment tensors of real data, can be produced from the P-to-S conversions at the interfaces between two different materials. Since these interfaces are usually unknown in the shallow part of volcano, single forces cannot be taken into account when interpreting results. On the other hand, De Barros et al. [2011] suggest that the single forces will accommodate most of the errors, so they should be included in the inversion procedure but should not be interpreted as part of the source solution. Including single forces in the solution appears to stabilize the results obtained using the 2011 data and therefore this paper focusses chiefly on results for which single forces are included in the inversion procedure. The single forces themselves are not interpreted as even though they can be physically realistic in volcanic environments, they appear contaminated by model errors.
Emphasis was also placed on the results which included single forces in the inversion procedure for the 2009 experiment [Eyre et al., 2013].

10.1002/2014JB011108
Constrained moment tensor inversion is also performed to test the robustness of the solutions, with the source mechanism constrained as (1) an isotropic source and (2) a crack with a grid search for the optimum dip and azimuth (for details, see Lokmer et al. [2007] and De Barros et al. [2011]). We reiterate that the single forces terms can accommodate model errors  and are therefore not physically interpreted. The number of data for every inversion is 3 components × 15 stations × N f (number of frequency points). The number of unknowns decreases from 9 × N f for the moment tensor + single forces inversion to 6 × N f (moment tensor only inversion) and 4 × N f for crack + single forces constrained inversion.
As moment tensor inversion takes into account the amplitudes recorded at each station, site amplification effects between stations must be removed to gain accurate solutions. These are measured using the RMS amplitudes of the coda waves of several teleseismic earthquakes that were recorded during each experiment on all stations [e.g., Zecevic et al., 2013], filtered between 0.3 and 1.3 Hz. The earthquakes used were located between 1000 and 1500 km from Turrialba Volcano and had magnitudes greater than 5.5 so that the coda for each event contained a measurable low-frequency component.
As the LP events do not appear to group well into families, moment tensor inversion is carried out for a significant proportion (107 events,~20%) of all of the LP events observed. The events chosen are those with optimal signal-to-noise ratios (generally the larger events), in order to give an accurate inversion solution.
We expect that inverting a large sample of LPs observed at Turrialba Volcano will yield a representative picture of all LP seismicity at Turrialba Volcano.

LP Location Results
Joint inversions and grid searches are implemented for 107 events, filtered between 0.3 and 1.3 Hz. An example of the results of the grid search location technique for a single event is shown in Figure 3. The figure shows three slices through the grid of possible source locations used for the inversions, with each grid point colored according to the residual obtained in the moment tensor inversion for that grid point. The vertical slices are taken through the minimum residual, while the horizontal slice is located at a depth of 2660 m asl (680 m below the summit). The most likely locations (the points with the lowest residuals) are shown in blue. It can be seen that the resolution in the horizontal plane is very good, with the location between the active Southwest and Central Craters, while the depth is less well constrained but is less than 800 m below the summit.   source area may be due to uncertainty in the medium velocity, but as these events do not belong to a family it is not surprising that the sources are scattered in a region. The results in the x-y plane are consistent with the results of array analysis and are also very similar to solutions obtained when single forces are not included in the inversion (Figure 4b), demonstrating the robustness of the solutions. However, the locations obtained when single forces are included in the inversion appear to cluster closer together, and there are a smaller number of outliers than for inversions excluding single forces, suggesting greater reliability in the results.

LP Source Mechanism Results
Moment tensor inversion solutions for the source locations with the lowest residuals are investigated by fixing the location and redoing the inversion. An example of the results from a single LP event when single forces are included in the inversion is shown in Figure 5. Note the pulse-like waveforms of the obtained moment tensor solution. The normalized waveform fits between the real and reconstructed data are also shown, showing good fits for the summit stations but poorer fits at more distal stations due to the intrinsic weighting of the stations by their recorded amplitudes. This is acceptable due to the strong oscillating path effects at stations farther from the source [Bean et al., 2014]. An eigenvector diagram is also shown, constructed using a method from Chouet et al. [2003] where the eigenvectors are calculated every 0.01 s when the amplitude of one of the moment tensor components is greater than 90% of the maximum amplitude. Finally, as described in section 3.1, singular value decomposition [Vasco, 1989] is used to extract a common source-time function with the scalar moment tensor containing its value for each component.
The eigenvectors are oriented with the long eigenvector 21°from vertical (this angle is herein named θ) and an azimuth anticlockwise from east (herein named φ) of 332°, with eigenvalues of 1.0:1.2:1.8. From the moment tensor decomposition of Vavryčuk [2001], this result is 74% isotropic. However, there is one larger eigenvalue and two shorter eigenvalues, suggesting that the mechanism can be described as a crack. This is demonstrated in the decomposition where compensated linear vector dipole (CLVD) components are found to contribute 18% to the solution. If the source is a pure tensile crack, then the long eigenvector is orthogonal to the crack plane; hence, θ gives the dip angle and φ the dip direction. The results, however, also show a double-couple component in both cases of 9%. It is therefore difficult to visualize the source mechanism. We reiterate that the single forces terms tend to "absorb" model errors , but the forces themselves are not reliable and are not constrained further.
In order to determine whether these complex results are unique to this one LP event, 107 events are inverted and summarized in Figure 6. The results of the moment tensor decomposition [Vavryčuk, 2001] (Figures 6a and 6b) and the orientations of the largest eigenvectors (θ and φ) (Figure 6c) are shown. Primarily, the figure shows results for inversions where single forces are included, but Figure 6a also shows results for moment tensor-only solutions for comparison and completeness. The results of the principal component analysis are shown as histograms in Figure 6a. The comparison between the results of inversions including and omitting single forces shows that similar distributions of results are resolved for each (i.e., similar peaks in the histograms), but the standard deviations in component percentages for inversions including single forces are much smaller. It is more likely that source mechanisms for LP events observed at the same volcano do not vary to a large degree, and therefore, these findings corroborate the work of De Barros et al.
[2013] who (as previously stated) suggest that moment tensor solutions from inversions including single forces are more reliable. This supports our decision to focus on results obtained when single forces are included in the inversions. Figure 6b shows the results including single forces as a triangular graph to show the full range of results. The triangular graph also shows the residual from the moment tensor inversion for each event using a color scale. A strong isotropic component is evident for all events, peaking at~75%. The CLVD and double-couple components are much less significant but do appear to contribute, with the CLVD component peaking at~5% but for some events contributing up to 40%, and the double couple peaking at~10% but for some events up to 45%. The variability here, together with the different source-time functions, is consistent with the LP events not presenting as families, i.e., events look different as the source mechanisms are different.
The mechanism orientations are shown in Figure 6c, both for the unconstrained inversions (MT + F) and for inversions when the source is constrained as a crack mechanism with a grid search on the orientation (Crack + F).   inversions. Although there is a large scatter in the results, it can be seen that for the majority of the events (especially when the residual is small) that θ is relatively small (30°or less). The cracks are therefore shallowly dipping, meaning that φ can be quite unstable. There also appears to be a preferential φ direction, with most events orientated with a φ of approximately 180°to 200°anticlockwise from east, i.e., due west, and perhaps a smaller cluster orientated roughly in the opposite direction (340°-360°, due east). Therefore, the cracks change dip around the horizontal position, from slightly eastward to the slightly westward dipping. Again, the variability in the orientations is consistent with the lack of similarity between event waveforms.

Joint Inversion and Grid Search Results: Discussion
The LPs are located at shallow depths below the volcano (<800 m). The results for the events with the lowest residuals, suggesting the most reliable inversions, mostly give depths between 700 and 800 m below the summit, although this may be several hundred meters shallower or deeper due to the uncertainty in the velocity model which can strongly affect the depth determination. However, the shallow results are in good agreement with LPs previously observed at Turrialba Volcano [Eyre et al., 2013] [Zecevic et al., 2013], and Kusatsu-Shirane, Japan [Kumagai et al., 2002;Nakano et al., 2003].
Most of the events appear to be closely clustered together below the active Southwest Crater. This is again in good agreement with results from 2009 [Eyre et al., 2013]. There are several events that appear to be outliers, however. It is most likely that these outliers are erroneous, as they generally have higher Journal of Geophysical Research: Solid Earth 10.1002/2014JB011108 residuals and are poorly constrained. The shallow locations below the active crater suggest a relationship with shallow processes and this needs to be taken into account when interpreting the mechanisms.
The moment tensor inversion results suggest that the source mechanisms are strongly isotropic. However, the mechanisms also appear to have small CLVD and double-couple components. As stated in section 3.1, the single forces components are likely to be unreliable and are not interpreted physically. These results suggest a very complex source mechanism. An isotropic and CLVD mechanism combined is interpreted as a crack mechanism. With a small contribution from double-couple included, the mechanisms could be interpreted as a crack mechanism with a small amount of shearing, akin to a transtensional or transpressional crack mechanism. A method for determining the mechanism parameters for this type of source is described in Vavryčuk [2001]. A tensile earthquake (an earthquake with tensile faulting or combining shear faulting and tensile faulting) can be described using a slip vector that is not restricted to orient within the fault plane and deviates from the fault plane, causing its opening or closing. This slip vector is labeled ⇀ u ½ and has an angle from the fault plane (labeled Σ) of α (Figure 7a).
The α and the λ/μ ratio (κ) at the source can be calculated using [Vavryčuk, 2001], where M Ã max and M Ã min are defined as the maximum and minimum deviatoric eigenvalues of the moment tensor solution and c ISO and c CLVD correspond to the components of isotropic and CLVD energy, respectively. The angle between the normal to the fault and the tension À ⇀ T Á axis is β = 45°À α/2, and this is also the angle between the slip vector and the tension axis. The angle (θ) between the tension axis and vertical is calculated using the eigenvectors of the moment tensor solution. From this, the orientation of the fault plane from vertical (δ) can be calculated using δ = θ -β.
The results are investigated for events with residuals obtained from the moment tensor inversion of less than 0.25 to remove events for which the mechanisms are more poorly resolved. Events with source locations that appear to be outliers are also removed. For the remaining 42 events the average of the double-couple components is 16%, of the isotropic components is 72% and of the CLVD components is 12%. The standard deviations are 8%, 8%, and 9 %, respectively. The relationship between α and the proportion of double-couple, isotropic, and CLVD components can be plotted, as shown here in Figure 7b [Vavryčuk, 2001] for the mean value of κ of 7.3. Low values of α correspond to a high degree of shearing as the slip vector is almost parallel to the fault plane, whereas high values of α correspond to tensile mechanisms. Using solutions with residuals lower than 0.25, the mean of the absolute values of α is 21°, while β is 34°, θ is 23°, and δ is À12°. The schematic geometry of the slip movement is shown in Figure 7a, with the relationships between α and the percentage of isotropic, double-couple, and CLVD components (Figure 7b). This shows that, even with a large value for the isotropic component, a small value of α can be found, i.e., the source mechanism does not deviate much from a shearing movement. Intuitively the results are very surprising as an α of 21°suggests a strong-shear component and a low tensile component, apparently at odds with the small double-couple component calculated from the decompositions of the solutions. Figure 7c shows a histogram of the distribution of α for all events inverted. There appears to be a large variability in α from À60°to 60°, and the standard deviation of α is 27°. Negative values of α suggest transpressional fracturing, and positive values suggest transtensional fracturing; however, due to potential errors from mismodeling the velocity structure, we do not believe that it is possible to accurately resolve the polarity of α. The low angles calculated for α therefore strongly support the possibility of a combined tensile-shear mechanism. Double-couple components have been observed in other studies of LP events [Cusano et al., 2008]. However, as the principal models of LP source mechanisms attribute them to fluid oscillation processes, these components are often interpreted as insignificant or explained as errors resulting from the calculation of the Green's functions. In such cases, it is possible to interpret the source mechanism as a crack, as was construed in the 2009 study of Eyre et al. [2013]. However, as the double-couple component is similar in size to the CLVD component, it would be difficult to justify including the CLVD component in the interpretation while ignoring the double couple. If the LP source is truly a tensile-shear mechanism, unidirectional slip motion should be observed in the source-time function, i.e., the band-limited response to a ramp function. On the other hand, a resonant source [Chouet, 1986] does not necessarily produce a source-time function with the low-frequency component of a ramp function. To test the hypothesis of whether the calculated source-time function can rule out a tensile shear mechanism, the source-time function obtained from singular value decomposition of the moment tensor solution for an LP event is plotted in Figure 8a, along with the band-limited signal (band-pass filtered between 0.5 and 1.3 Hz) produced by the ramp function shown in Figure 8b which has a rise time of 0.5 s. As expected, the onset and overall shape of both signals are similar; however, the LP source-time function shows some extra higher frequency fluctuations later in the signal. These could be caused by errors in the inversions such as mismodeling of the velocity structure, leading to path effects mapping into the source-time function. Previous near-field observations have shown that the codas of recorded seismic events in volcanic environments are strongly influenced by propagation effects [Bean et al., 2008[Bean et al., , 2014Harrington and Brodsky, 2007;Lokmer et al., 2007] which will be present in the obtained sourcetime functions. Another issue is that details of the source-time function are not possible to capture for such small events, where the amplitude of the very long period part (which would be produced by a static shift) is below microseism level, the shallow path effects are strong, and the seismic phases are intertangled. Furthermore, as the mechanism has a tensile component, the source-time function is likely to be more complicated than in the case of simple shear. The opening of a crack is followed by a closing, as the pressure within strongly decreases during opening [Eaton et al., 2014]. This means that tensile mechanisms only exhibit a small unidirectional movement. In the case of LP events at Turrialba Volcano, as the source is complex, involving both the opening and closing plus shearing of the crack, the source-time function is also more complex than for the simple shearing case, which may explain the source-time function observed. Due to all of the factors discussed, we can only consider a match by comparing the broad-scale shape of the source-time function (which includes the static component) to the filtered ramp function, and especially the beginning of the signal as the coda part is the most likely to be affected by the factors discussed. Therefore, although the observed source-time function does not fully match that of a unidirectional source, it is similar enough to show that a tensile-shear mechanism cannot be ruled out.
The obtained κ values should provide an indication of the fractured material. The mean κ value of 7.3 corresponds to a Poisson's ratio of~0.44. Such high values suggest that the source region either consists of (1) highly fractured or porous media saturated with fluids (as P wave velocities are driven by the fluids, while S wave velocities depend on the solid frame) or (2) ductile media (which may or may not be caused by high temperatures). For example, dense sands (which have a ductile response to stress and are likely highly similar to the unconsolidated material typically found at the near surface in the summit regions of volcanoes) have a Poisson's ratio of 0.45 [Greaves et al., 2011]. As fluids, structural heterogeneities, unconsolidated materials, and high temperatures are likely to be present in the shallow part of the volcano, Figure 8. (a) Comparison between the source-time function obtained from singular value decomposition of the moment tensor solution of an LP event from Turrialba Volcano, filtered between 0.5 and 1.3 Hz (blue) and (b) the signal obtained by filtering the ramp function with a rise time of 0.5 s in the same frequency band (red). Note the similarities between the two signals: the first peak and overall shape. Higher-frequency fluctuations observed later in the decomposed signal may be caused by errors in the inversion procedure or may be due to the complex interaction within the source of unidirectional shearing plus tensile opening and closing.
Journal of Geophysical Research: Solid Earth 10.1002/2014JB011108 either of these are strong possibilities. However, it must also be taken into account that the obtained value of κ can be strongly influenced by the accuracy of the velocity model.

Scaling Laws of LP Events
The results outlined above show that the events are caused by a tensile mechanism with a likely strong-shear component. These results can be well explained by the new source model for LP events of Bean et al. [2014] where the events are generated by slow failure in the shallow edifice. However, further analysis of the LP events is necessary to investigate how well the model corroborates the data and to explore brittle failure as a valid model. Tectonic earthquake seismology uses empirical scaling laws to link the seismological observations to the source properties, and the above observations encourage us to apply such analysis to these LP data to further investigate their source.

Frequency-Magnitude Relationship
The frequency-magnitude relationship [Gutenberg and Richter, 1944;Ishimoto and Iida, 1939] usually gives a b value of approximately b = 1 for tectonic earthquakes. Figure 9 shows a plot of the LP event maximum amplitudes (plotted logarithmically to use as a proxy for magnitude) against the cumulative number of events for the LPs observed at Turrialba Volcano during the 2011 experiment, recorded on station TCCR (which has a higher signal-to-noise ratio than WCR1). Amplitude can be used as a proxy for magnitude only for earthquakes which are at the same distance from the station, which is satisfied here. It is noted that the slope of cumulative number versus amplitude in log-log plot is actually known as the m value [Ishimoto and Iida, 1939]. Above the visually estimated detection threshold of~3 × 10 À6 m/s the relationship appears to show nonpower law scaling with a deficit of larger events.
One possible interpretation for this behavior is that the steepest gradient (À2.9, line shown in Figure 9) corresponds to the true b value, and the real detection threshold is~2 × 10 À5 m/s. Note that here we assume that amplitude is proportional to the magnitude (b = m), but even if we assume a relationship of m = b + 1 instead [Suzuki, 1959], we obtain b = 1.9 which is still higher than that observed for ordinary earthquakes. The smaller events (amplitude < 2 × 10 À5 m/s) do not fit the relationship due to the variable levels of noise (e.g. volcanic tremor) resulting in an incomplete catalogue. Other potential issues include the relatively small size of the catalogue and the small magnitude range of the catalogue (~2 orders of magnitude). For amplitudes below~3 × 10 À6 m/s, the graph is flat as the events are no longer detectable. In this interpretation larger LPs occur less often than expected based on the typical earthquake scaling law. However, large b values have been observed for earthquakes and volcano-tectonic (VT) events in volcanic areas, e.g., at Mount St. Helens and Mount Spurr [Wiemer and McNutt, 1997], Mount Pinatubo [Sánchez et al., 2004], and the off-Ito volcanic region [Wyss et al., 1997], so this is also feasible for LPs. Warren and Latham [1970] state that b values up to 3 are possible in earthquakes induced by high thermal gradients, supported by Wiemer and McNutt [1997]. Another possibility is that the high b value could be caused by high heterogeneity of the material where the fracturing occurs [Mogi, 1962] and/or a low stress regime [Schorlemmer et al., 2005], both of which are highly likely at shallow depths on volcanoes. Finally, Wiemer and McNutt [1997] also propose that increased pore pressure due to heating of the groundwater could increase pore pressure and so lower effective stress, which is also very likely at shallow depths beneath the Southwest Crater at Turrialba Volcano. The high b value of the LPs could therefore be caused by a combination of these factors.
Nonetheless, it is felt that the section of the graph with a lower gradient occurs over a much larger range of amplitudes than would be expected if it was caused by incompleteness of the catalogue. Therefore, the behavior is interpreted as a real manifestation caused by the source processes of the events. Bean et al. [2014] observed a similar nonpower law scaling for failure at the brittle-to-ductile transition in materials similar to those observed at the near surface of volcanoes simulated using a damage mechanics model. . LP event amplitudes (plotted in logarithmic scale as a proxy for magnitude) against cumulative number of events greater than that amplitude. The catalogue appears to be complete above the visually estimated detection threshold of~3 × 10 À6 m/s. The line of best fit for the steepest section of the plot has a gradient of~À2.9.
Journal of Geophysical Research: Solid Earth 10.1002/2014JB011108 They also stated that this relationship is very similar to magnitude-frequency relationships of LPs observed at Mount Etna, Sicily, and at Turrialba Volcano during the 2009 experiment. The observations made here are consistent with these earlier findings.

Corner Frequency Analysis
Earthquakes follow a corner frequency (f c ) versus magnitude (M 0 ) scaling. The most commonly used relationship is M 0 ∝ f c À 3 [Aki, 1967;Geller, 1976;Hiramatsu et al., 2002;Kanamori and Rivera, 2004]. Before this analysis is performed, we must first confirm that the LP spectra appear similar to earthquake spectra, with flat spectra up to a corner frequency where the energy begins to fall off (the fall off is often modeled using an ω À2 model). The unfiltered displacement spectrum for the event from Figure 2a is plotted in log-log scale in Figure 10. This shows that the events do behave like typical earthquakes, and a corner frequency can be picked. As previously mentioned, the peak above 10 Hz corresponds to continuous noise (spectral estimates for pre-LP signal noise are plotted in blue).

Synthetic Tests
A concern with implementing corner frequency analysis on LP events at volcanoes is that the LPs are recorded in the near field and have a tensile component, since the corner frequency concept is valid in the far-field approximation for double-couple mechanisms. On a theoretical basis, Eaton et al. [2014] show that the corner frequency relationship is still valid for events with tensile components. To investigate whether the corner frequencies can be accurately recovered in the near field, synthetic tests are conducted.
The Stokes solution of elastic wave radiation [Aki and Richards, 2002;Lokmer and Bean, 2010] is used to kinematically calculate synthetic seismograms. The near-field and intermediate-field effects can be included in the calculations. In order to simulate finite fault ruptures, the source comprises a 2-D grid of "tensile crack" point sources representing the rupture plane. These point sources rupture unilaterally along the longer dimension of the rupture plane. In order to yield a range of event magnitudes, the rupture plane ranges in size from 50 m to 1000 m long (with the width equal to half of the length).
To test whether the relationship should hold specifically for the LP events recorded at Turrialba Volcano with the stations used, the rupture plane is located in a similar location to the recorded LP events (UTM 196,920 m east, 1,108,760 m north and 2660 m asl) and has a similar orientation with θ = 20°and φ = 190°. The first station tested is WCR1, which has the highest amplitudes and least propagation effects for the events as it is closest to the source and therefore would be the optimum station to employ for source property studies.
In order to find the corner frequencies, the spectra of the synthetic events are fitted to the ω À2 model using the equation: where f is the frequency and f c is the corner frequency [e.g., Abercrombie, 1995]. The fit is implemented for a range of corner frequencies to find the corner frequency with the optimum least squares fit between the equation and the data. The fit also outputs the amplitude for the flat part of the spectrum of the event, which is used as a proxy for the magnitude.
The relationship between corner frequencies and magnitudes for the synthetic events in the far field is shown in Figure 11a, along with the fits between the spectra and fitting function (Figure 11b). A strong inverse Figure 10. Displacement frequency spectrum of the LP event from Figure 2a recorded at WCR1, in log-log scale (black). The noise spectrum taken before the event is also shown (blue), demonstrating that the event spectrum is above the noise. The shape of the spectrum appears similar to that of a regular earthquake, allowing a corner frequency to be picked. The peak above 8 Hz corresponds to continuous noise and so is ignored in the corner frequency analysis.
Journal of Geophysical Research: Solid Earth 10.1002/2014JB011108 relationship between the two is demonstrated as expected. The relationship is not quite f c À3 as is expected (the slope of the best fit line is À3.45 instead of À3). When the near-field (and intermediate-field) term is included, invalidating the far-field assumption, the results are worse, but an inverse relationship can still be recognized with a slope of À5.03 (Figure 11c). To check that this is a real problem caused by the close proximity to the source, the analysis is also carried out at a station further from the source, station DIVI. Here the analysis results in gradients of À3.13 and À2.57 excluding and including the near-field term, respectively, much closer to the expected solution of À3. Figure 11. (a and c) Relationship between corner frequency and magnitude for synthetically created events located in a similar location to the real events and recorded at station WCR1. Also shown are the respective fits between (b) the synthetic event spectra and (d) the function used to estimate the corner frequencies.
L denotes the length of the rupture plane. The results for the far-field term only in Figures 11a and 11b and the results including the near-field term in Figures 11c  and 11d. The best fit lines are shown: slope of À3.45 in Figure 11a and slope of À5.03 in Figure 11c. For ordinary earthquakes, À3 is expected. The results suggest that the corner frequencies of events are affected by recording close to the source, especially by the near-field term; however, an inverse linear relationship is always present.
Journal of Geophysical Research: Solid Earth 10.1002/2014JB011108

Analysis of LP Data
The synthetic analysis has shown that ideally a far-field station should be used for the analysis in order to reduce the near-field effect. Unfortunately, this is not possible as these stations have much larger propagation effects included in the LP spectra [Bean et al., 2008] which can strongly affect the results of seismic source studies [Giampiccolo et al., 2007] have a lower signal-to-noise ratio and have a smaller range of LP event amplitudes. Therefore, station WCR1, on the rim of the active Southwest Crater, is used for the analysis. The synthetic tests have demonstrated that the gradient of the best fit of the data cannot be interpreted, but if an inverse relationship is uncovered, this still has implications for the source mechanism.
As the LPs recorded are relatively small (of the order of 1 × 10 À5 m) in comparison to the earthquakes usually analyzed using this method and the noise is quite high (as shown in Figure 10), it is very difficult to determine the corner frequencies of the events. In order to overcome this issue, the following procedure is implemented. The events are first extracted and cut into 10 s long traces before integrating into displacement. Their means and trends are removed, the traces are tapered using a Hann window and are then transformed into the frequency domain. The spectra are smoothed using a moving average. Corner frequencies of the events can now be obtained using the same technique used for the synthetic data. Fits are implemented for a frequency range of 0.3 to 7.0 Hz, just above the microseismic noise range.
The results of the analysis are shown in Figure 12a. Figure 12b shows examples of the fit of the ω À2 model used to calculate the corner frequencies with the data spectra of random individual events. It can be seen that there does appear to be a relationship between the corner frequency and the magnitude of the events. The L1 norm best fit has a gradient of À2.0 rather than the À3.0 expected in the dry, brittle fracture of earthquakes when analyzed in the far field. However, as the noise and the distribution of the event sizes can bias the linear regression, we computed two L1 norm linear regressions: the corner frequency as a function of the amplitude and the inverse relationship. The average value we obtain for the gradient is À2.9, which is very close to the À3.0 expected. Therefore, the f c À3 relationship is consistent with the data. The observed inverse relationship is consistent with a brittle failure origin for the Turrialba Volcano LP events.

Discussion and Conclusion
The observations and results presented in this paper corroborate the possibility that LP events at Turrialba Volcano may be generated by a failure mechanism. First, LPs observed close to the source (i.e., at the summit stations) are often pulse like, and this is evidence that especially conflicts with the Chouet [1986] model of fluid-filled cavity resonance. This fact could be explained by a high damping of the resonance due to a high-viscosity fluid or low-crack wall stiffness [Chouet, 1986], but it seems suspicious that a model to describe resonance should apply for pulse-like waveforms. On the other hand, a material failure model does not require such specific tuning.
Second, the frequency-magnitude relationships and corner frequency-magnitude relationships are very similar to those observed in rock fracture experiments and earthquake studies. More work is required to analyze whether other possible source mechanisms or triggers can cause the relationships observed, but this goes beyond the scope of this paper. Together though, these observations appear to suggest that a mechanical failure mechanism is more likely than a resonant mechanism for LP event generation at Turrialba Volcano.
Finally, the results from the moment tensor inversion suggest a tensile mechanism with a strong-shear component. This is difficult to explain using the Chouet [1986] model but can readily be explained using the model of Bean et al. [2014] as tensile (mode I) cracks combined with shearing, which results in transtensional-or transpressional-like faulting. Transtensional faults seem especially plausible in volcanic environments due to flank instabilities and hydrofracture in hydrothermal systems, and transpressional faults could be caused by gravitational collapse. Dipping slump faults have been mapped at volcanoes, for example, Stromboli [Apuani et al., 2005], and shallow-dipping detachment faults have been proposed at Mount Etna [Rust et al., 2005]. In the case of Stromboli, a VLP event was constrained as a slow-slip double-couple event caused by dip slip [Cesca et al., 2007]. This type of faulting fits well with the orientations of the mechanisms obtained for the LPs at Turrialba Volcano, as they have shallow dip angles. Turrialba Volcano sits within a SW-NE trending graben structure with a strike-slip component of motion [Soto, 1988]; however, linking the tectonic regime to the LP mechanisms is very difficult in the shallow part of the volcano as the stress regime can be completely altered by the presence of fluids, topography effects, gravitational instability, and structural heterogeneities near the summit. Figure 12. Relationship between corner frequency and magnitude for LP events recorded at Turrialba Volcano for station WCR1. (a) The results from the corner frequency analysis of the LP events. As the magnitude range is small and the fits are not exact, the events are quite clustered, but there does appear to be a clear inverse relationship between corner frequency and magnitude. A line of best fit (calculated using the L1 norm) is plotted through the data (red line), with a slope of À2.04. A second line of best fit has also been calculated using the inverse of the relationship, as the noise and the distribution of the event sizes can bias the linear regression, and is plotted as a dashed red line, with a gradient of À3.76. The theoretical À3 earthquake relationship is also plotted (black with markers) and is also consistent with the data. (b) Examples of the fits between the smoothed spectra of nine random LP events (blue) and the fitting functions (red).
The fits appear to be reasonable.

10.1002/2014JB011108
A likely underlying cause of the events observed at Turrialba Volcano therefore appears to be deformation. The majority of the events are located below the active Southwest and Central Craters, as was also constrained for events recorded in 2009 [Eyre et al., 2013]. This suggests that small-scale deformation is occurring below these active craters. However, the recovered source-time function is not accurate enough to distinguish between inflation, which could be caused by an influx/rising of magma or hydrothermal processes, and deflation, which could be due to the loss of large quantities of gas and edifice collapse. As previously mentioned, large-scale deformation has not been observed at Turrialba Volcano since 2007 [Martini et al., 2010] so the deformation causing LP events is likely to be small and localized beneath the active craters.
Future work is necessary to explore the slow-failure model for LP event generation as described by Bean et al. [2014]; however, the observations and analysis described in this paper suggest that brittle failure is a possible LP generating mechanism on Turrialba Volcano.