Monitoring methane emission of mud volcanoes by seismic tremor measurements: a pilot study

A new approach for estimating methane emission at mud volcanoes is here proposed based on measurements of the seismic tremor on their surface. Data obtained at the Dashgil mud volcano in Azerbaijan reveal the presence of energy bursts characterized by well-determined features (i.e. waveforms, spectra and polarization properties) that can be associated with bubbling at depth. Counting such events provides a possible tool for monitoring gas production in the reservoir, thus minimizing logistic troubles and representing a cheap and effective alternative to more complex approaches. Specifically, we model the energy bursts as the effect of resonant gas bubbles at depth. This modelling allows to estimate the dimension of the bubbles and, consequently, the gas outflow from the main conduit in the assumption that all emissions from depth occur by bubble uprising. The application of this model to seismic events detected at the Dashgil mud volcano during three sessions of measurements carried out in 2006 and 2007 provides gas flux estimates that are in line with those provided by independent measurements at the same structure. This encouraging result suggests that the one here proposed could be considered a new promising, cheap and easy to apply tool for gas flux measurements in bubbling gas seepage areas.


Introduction
Mud volcanism and diapirism are well-known phenomena occurring whereby fluid-rich fine-grained sediments ascend within a lithologic succession by buoyancy due to the bulk density contrast between an overpressured muddy mass and an overburden one of higher density.Such overpressure is also the effect of gas production at depth (mainly methane) induced by organogenic activity within water saturated sediments buried during subduction or rapid subsidence phenomena (see reviews by Higgins and Saunders, 1974;Kopf, 2002).These structures have been extensively studied in the last years because they could be a marker of hydrocarbon reservoirs and for hazard induced by violent mud extrusion during drilling (Bagirov et al., 1996;Kopf, 1999;Robertson et al., 1996;Davies et al., 2010;Rudolph et al., 2011).
Interest in monitoring mud volcanoes activity and in particular in their gaseous emissions comes from recent data, which suggest that methane emission at these structures significantly contributes to the global budget of greenhouses gases (Etiope and Ciccioli, 2009;Judd, 2005).Such monitoring is also of interest for hazard assessment in many countries where dramatic inland mud eruptions have been historically documented (see, e.g.Aliyev et al., 2009;Mazzini et al., 2007) and as a proxy for co-seismic strain detection (Albarello, 2005;Kopf et al., 2010;Martinelli and Panahi, 2005).
In principle, monitoring the activity of mud volcanoes in terms of gas outflow is possible by capturing emissions at the surface of emitting conduits (Kopf et al., 2010;Martinelli et al., 1995).However, as this approach works well when small structures are of concern, when several independent seepage structures and diffuse microseepage exist (e.g.Etiope et al., 2011), the overall gas outflow can be only indirectly estimated in this way.In addition, significant logistic difficulties Published by Copernicus Publications on behalf of the European Geosciences Union.also exist which make this kind of measurements infeasible and expensive in many contexts (Kopf et al., 2010).These difficulties explain also why most estimates of gas outflow from inland mud volcanoes are only performed by spot measurements (e.g.Etiope et al., 2004a, b;Etiope, 2005) A possible way to face these problems is the use of indirect estimates provided by seismic ground motion observations, which revealed to be of great importance in the monitoring of "hot" volcanoes (e.g.Chouet, 2003).These observations indicate that, in the case of hot volcanism, specific seismic signatures are present that are strictly linked to magma ascent mechanism and are the result of a complex interaction between magma fluids and surrounding bedrock (e.g.Konstantinou and Schlindwein, 2002).Several dynamic models were proposed to link seismicity and physical phenomena occurring in the magma reservoir and in the conduit (e.g.De Lauro et al., 2008 and references therein).In particular, energy bursts commonly observed in the infrasonic range (few hertz to tens of hertz ), both in acoustic and seismic measurements, were associated with bubble dynamics (resonance or breaking) inside the volcano reservoir (e.g.Ripepe and Gordeev, 1999;Vergniolle and Brandeis, 1994;Vergniolle et al., 1996;Vergniolle and Caplan-Aurebach, 2004;De Martino et al., 2011, 2012).
It is well known that bubbling also plays an important role in mud volcanism.Gas bubbles have been recognized also in soft sediments (e.g.Best et al., 2004), and they are important in the study of gas hydrate dynamics (e.g.Bourdeau et al., 2005).Low permeability of clays in mud-volcano areas (Kopf, 2002) suggests that, in the lack of large mud outflow (typical of quiescent phases), gas propagation from the reservoir should mainly occur by the uprise of bubbles (Etiope and Martinelli, 2002;Albarello, 2005).Actually, visual monitoring of bubbling at the surface of emitting vents was considered for estimating gas outflow at mud volcanoes (Etiope et al., 2004a).The basic limitation of this last approach is that it is feasible when strictly localized emission occurs, without any other dry seepage.A different approach to the problem of counting of bubbles rising from submerged vents is based on acoustic emissions induced by active hydroacoustic echosounding (Greinert and Nützel, 2004;Greinert et al., 2006;Nikolowska and Schanze, 2007;Nikolowska et al., 2008).Natural acoustic emissions (Leighton, 1994) from bubbles at depth can be also considered for monitoring inland structures with the hypothesis that they provide seismic signal that can be detected at the surface, as in the case of hot volcanoes.In principle, this approach could allow better estimates of the overall degassing activity of volcanoes, since outflow is detected at depth before the partitioning of gaseous emissions at different surface seepage structures.
In this framework, seismic measurements at the Dashgil mud volcano in Azerbaijan were carried out in the frame of the NATO-CLG program "Remote Sensing to Monitor Methane Emissions in Mud Volcanoes" (Ref.982053) during two field surveys in 2006 and 2007 respectively.Dashgil volcano has been extensively studied in the past, and thus it represents an ideal case study to test new measuring tools and procedures.In the following, major features of Dashgil volcano are described at first.Then, seismic data processing is outlined along with major results and possible interpretations.

The Dashgil mud volcano
Dashgil mud volcano is a 2-km-wide dome located in Azerbaijan near the western coast of the Caspian Sea (Fig. 1).The bulk of Dashgil is formed by clayed mud with variable amounts of clasts in the matrix originating from the Maikop formation at several kilometres depth (Kopf et al., 2010).Dashgil volcano is one of the most active mud domes of the Earth, and the last major eruptions occurred in the years 1882, 1902, 1908, 1926and 1958(e.g. Jakubov et al., 1971;;Guliyev, 1992;Aliyev et al., 2009).Beyond this paroxystic activity, it also represents a good example of mud volcano with high seep activity in the dormant period (Planke et al., 2003).Due to these features, Dashgil volcano has been extensively studied in the past from the geomorphologicalgeostructural (Hovland et al., 1997;Roberts et al., 2011), geochemical (Planke et al., 2003), and geotechnical (Kopf et al., 2009) point of view and was the object of long-term monitoring of gaseous emissions (mainly Methane) at the summit area (Kopf et al., 2010).
The summit area of the Dashgil volcano (Fig. 1) presents the typical structure of a mud volcano summit caldera (Evans et al., 2008).Gas, mud and water seepage occur at a number of different structures (Hovland et al., 1997;Mazzini et al., 2009), and their distribution is probably controlled by structural features (Roberts et al., 2011).In the southwestern sector of the summit area, a north-south oriented field of active gryphons and mud cones exist.These last ones consist of few meter high trunk-cone structures filled or partially filled of bubbling clay.Gryphons instead are smaller seepage holes where no mud seepage occurs.At these structures, there is occasionally a deep bubbling sound and a faint oil smell inside at least some of the cones.On the basis of geomorphological considerations, Roberts et al. (2011) suggest that circular rim surrounding the gryphon area nearly corresponds to the main vent of the volcano.
Towards the eastern edge of the caldera, two lakes of low viscosity mud (salsas) exist.In both, gas bubbles rise from the bottom of the lakes, which reach maximum depths of about 10 m and 8 m, respectively (Kopf et al., 2009).Geochemical analyses (Planke et al., 2003) suggest that the water-dominated salsas represent expulsions of shallow waters, whereas deep-seated waters percolate slowly to the surface and are expelled as mud-water mixtures in gryphons and springs.Mazzini et al. (2009) suggest that a mixing of meteoric water with a deep originated brine could justify some peculiar characteristics (low Boron content) of one of the Dashgil's salsas (Salsa A).However, the studied area is in the BSh Köppen-Geiger climate area (Peel et al., 2007) and a water surplus able to recharge a phreatic aquifer is almost impossible.Furthermore, no gravel or sandy sediments able to help water infiltration were found in the Dashgil's area.During a field survey carried out in the frame of the NATO-CLG project (Ref. 982053) in 2007, two of the Authors of the present paper (M. and A.) found 18 O/ 16 O and D / H of +6.24 ‰ and −33.8 ‰ in the Dashgil gryphon and +3.09 ‰ and −30.9 ‰, respectively in the salsa (unpublished results), which were interpreted as the result of sediment exchange with sea-water-originated fluids and sediments confirming the connate character of the emitted fluids in the liquid phase.In this view, overpressured mud, water oil and gas rise up in a nearly vertical pipe from a depth of about 2-3 km (Jakubov et al., 1971;Hovland et al., 1997). Similarly, Hovland et al. (1997) suggested that all the structures in the crater are expression of a unique reservoir at depth.This suggestion was also supported by direct observations provided by Kopf et al. (2010).
Several attempts have been provided to estimate gas seepage at the Dashgil volcano.Earlier estimates were provided by Jakubov et al. (1971) and Hovland et al. (1997) suggesting flow rates of the order of 2×10 3 m 3 h −1 and 0.09 m 3 h −1 , respectively.Guliyev and Feizullayev (1997) provide estimates of emission rates in the wide range of 8 m 3 h −1 to 2 × 10 3 m 3 h −1 .Etiope et al. (2004b), by taking into account both venting (pools and gryphons) and microseepage, estimated a total average methane output around 140 m 3 h −1 .An estimate of the order of 1 m 3 h −1 was provided by Kopf et al. (2009Kopf et al. ( , 2010) ) after four years of nearly continuous gas monitoring at the salsas, i.e. a selected single emissions among the ones present in area.These measurements showed that gas flux is characterized by valve-type behaviour and episodically violent degassing: this could explain the large differences provided by spot-like measurements.
At the bottom of salsas, Kopf et al. (2009) identified the presence of conduits where overpressurized fluids outflow with an excess pore pressure very similar to the estimated lithostatic pressure.Morphological indications of possible updoming of the crater in association with gas outflow variations (Kopf et al., 2009(Kopf et al., , 2010) ) suggest that a shallow fluid reservoir could exist under the crater, where dynamic equilibrium exists between pressure of uprising gas and overlying dry clays.Gas present in this reservoir could be the source of both gas emissions at the salsas and of high viscosity mud outflow and bubbling at the mud cones.

Seismic tremor surveys
In order to minimize logistic troubles, the monitoring of mud volcanic tremor was carried out by using a portable seismometer, i.e. "Tromino".It is an ultra-compact (1 dm 3 ) and ultra-light (< 1 kg) instrument for high-resolution digital seismic noise measurement (see www.tromino.itfor details).Tromino is an all-in-one (no external cables exist) tri-directional, high-sensitivity seismograph equipped with a GPS apparatus for timing and location.The ground coupling is obtained through three spikes fixed in the ground.Its very low energy consumption allows long-lasting measurements (days) with internal supply (2 AA 1.5 V batteries) only.The data storage is internal (up to 2 GB).
Three sessions of seismic measurements were carried out in the crater area.The first two sessions (hereafter I and II) were performed on 23 June 2006 around noon (local time) near the mud cone and pool respectively (Fig. 1).Each session lasted 20 min and a sampling rate of 128 Hz was adopted.During measurements, wind was nearly absent and no visible source of anthropic noise was detected.Despite the fact that measurements were carried out few tens of minutes one from the other and few hundreds of meters apart, the two time series appear to be quite different to visual inspection (Fig. 2).Ground motion amplitude at the mud cone is about one order of magnitude larger than at the salsa.This is the effect of energy bursts (hereafter "transient events" or TEs) that are present in the first series and absent in the second series, even if it also concerns the background tremor.The shape of TE (see the inset aa in Fig. 2) strongly mimics the shape of bubble bursting events detected from acoustic measurements at Stromboli volcano in Italy (e.g.Ripepe et al., 1996;Vergniolle et al., 1996).
These encouraging outcomes suggested to repeat seismic measurements in a subsequent survey one year later.On the 1 July 2007, a new measurement session (hereafter Series III) was performed near mud cones (see Fig. 1).In this case the session lasted 20 h, during which the tromograph was left in a small hole drilled in the mud (20-30 cm of depth).Seismic tremor measurements started at about noon (local time) with a sampling rate of 128 samples per second.In the first half of the measurements session, a strong wind (typical of the area under study) affected the area providing a strong disturbance for most of the time.Despite of this, TEs were clearly observed again (see Fig. 2cc) confirming the possible presence of phenomena reflecting persistent activity in the summit area of mud volcano.
These findings suggested to focus our analysis on the TEs of the series I and III.In particular, due to the low level of disturbances of the series I, data processing procedure was set up on this series and applied subsequently on the series III.

Series I (June 2006)
As pre-processing, the series was frequency filtered with a high-pass filter with corner frequency equal to 1 Hz, in order to reduce the effects of the microseismic noise induced by sea wave activity (the Caspian Sea is few km aside).Then TEs were focused on.All the detectable TEs were picked by an automatic process in which we evaluate the maximum of the signal amplitude within non-overlapping windows of 1.2 s and compare the maxima of two adjoining windows.In this way an energy burst is identified when both the following conditions are satisfied: (1) the ratio between the maximum of the latter window and that of the former window exceeds a threshold equal to 1.4, and (2) the ratio between the maximum of the latter window and the standard deviation of all the signal exceeds 2. The optimal window length and the thresholds were empirically determined.The detected TEs were then visually checked looking at their waveforms to exclude the picking of spurious signals.Moreover, we calculated also the spectrogram of the signal in order to highlight possible low-amplitude TEs that could not be detected by the automatic picking process (less than ∼ 2 % of the events).In this way we detected 136 TEs during the 20 min of registration.
The selected TEs show simple waveforms constituted by a few cycles of a nearly monochromatic signal (see Fig. 2aa).In order to investigate the stability of the source process, we perform the spectral analysis and the cross-correlation analysis on the TEs.Specifically, we compute the fast Fourier transform (FFT) of each TE.In addition, we estimate the maximum of the cross-correlation function between one master TE and all the others, and several possible master TEs are used for the analysis, all of them providing the same results.
In this way two classes of TEs (Class A and Class B hereafter) were identified on the basis of the cross-correlation analysis, which provides a bi-modal distribution of its maxima independently of the master TE chosen.Moreover, these classes display TEs with different spectral content, although both are constituted by nearly-monochromatic events.In Fig. 3 the stacked waveforms of the TEs for each class are reported along with the corresponding mean power spectrum.TEs of the Class A show sharp spectral peak at about 3 Hz, whereas for Class B TEs the peak is in the range 4.5-5 Hz.Finally, 60 events of the Class A and 76 events of the Class B were detected.
In order to infer constraints on the source positions and on the source mechanisms driving the TEs, a polarization analysis (Kanasewich, 1981) was performed.In this procedure diagonalization of the covariance matrix derived from the three ground motion components (vertical, north-south, east-west) was achieved.The eigenvector corresponding to the highest eigenvalue is considered as the best estimate of the polarization vector.This vector is defined by two angles: the azimuth and the dip.The former measures the angle with respect to the reference direction of the seismometer (see Fig. 1); clockwise, the latter is the angle formed with the vertical direction.Another parameter called rectilinearity (RL) defines the degree of elongation of the oscillations: it is equal to 1 when polarization is perfectly linear, and it is 0 when the eigenvalues of the covariance matrix are equal and a preferential direction of oscillation cannot be defined.Only solutions with high RL (> 0.6) can provide reliable value of azimuth and dip angles, as in that case the oscillations are elongated enough to define a dominant polarization direction.
In Fig. 4 the time evolution of the three polarization parameters for the two classes of TEs is plotted.Here, the solutions of all the TEs of each class are stacked on.This visualization highlights the common behaviour of the polarization parameters for all the TEs and enforces the distinction of these classes of TEs.Moreover, it emphasizes the occurrence of the TEs within the background signal.In detail, we note very high values (> 0.8) of RL for both the classes, suggesting very elongated particle motions in correspondence with the TEs.In particular, RL of Class A events is higher than that of Class B on average.The two classes share stable dip centred at about 60 • (from the vertical), whereas azimuth angles are equal respectively to about 10 • and 30 • (from the reference direction in Fig. 1).Main features of TEs in series I are summarized in Table 1.

Series III (July 2007)
Starting from the results obtained for series I, a similar analysis on the series III was attempted.This is much longer than Table 1.Summary of TEs' main features in the two time series of tremor measurements at the Gryphon area.For each Class and sub-class of events, characteristic frequency (f ), polarization dip ( • from the vertical), polarization azimuth ( • from the reference direction), number of TEs in the relevant series (N ), time rate of TEs (number per hour) and radius estimate of the bubble responsible for the relevant TE are reported.series I but is affected by strong disturbances mainly due to the intense wind.Due to such a strong background noise, it is sometimes difficult to individuate the TEs on this series, above all for the low-amplitude TEs.To enhance TEs signal, which are characterized by high rectilinearity, a polarization filter was applied to the signal before the TEs' detection process.In particular, a 1.2 s sliding time window scanned the whole series.In each window, the RL parameter was computed and used to modulate the recording (see Jurkevics, 1988;Vidale, 1986).On purpose, modulation was performed by using a factor RL 4 .An example of the final series is plotted in Fig. 2c.After this pre-processing, the picking of the TEs was achieved as described above for the series I.In this way, about 1300 TEs were selected.Then, the FFT of each TE was computed to estimate the highest spectral peak.The distribution of these spectral peaks is plotted in Fig. 5 and shows the clear presence of three relative maxima.
In this case, differently from the series I, TEs with spectral peak close to 3 Hz are absent.Moreover, the distribution displays three clear peaks at about 4.5 Hz, 5.5 Hz, 6.5 Hz, thus  suggesting the existence of a sort of sub-classes of class B TEs. Thus, we separated the TEs in three kinds of TE on the basis of their spectral peak: 4-5 Hz, 5-6 Hz, and 6-7 Hz, hereafter labelled as sub-classes a, b and c, respectively.Examples of the waveforms relative to the three sub-classes are plotted in Fig. 6.
As a whole, 1120 events of this kind were detected (about 82 % of the whole set of "events" actually detected), in particular, 532 TEs of sub-class a, 380 of sub-class b and 208 of sub-class c.In performing cross-correlation analysis, the above subclasses were accounted for.Hence, we estimated the maxima of the cross-correlation function between a master TE chosen among the sub-classes a, b, c and all the TEs.This analysis was also repeated by considering as a master event the waveforms detected in series I both considering class A and class B. In Fig. 7a-c  The distributions in Fig. 7a-c display one maximum close to 0.7-0.8 and a long tail extending up to values lower than 0.5.This pattern confirms the existence of a variable source mechanism, which is not limited to the existence of one or more well-distinguished classes of events.In this sense, the distinction of sub-classes based on the spectral content of similar waveforms appears to be a good choice.The same behaviour is detected also for Fig. 7e, indicating that TEs of series I are contained in the series III, although this latter case displays a broader variability in the source process.On the other hand, the distribution of Fig. 7d shows one maximum at about 0.2-0.3 which suggests that TEs of class A do not exist in series III.
A polarization analysis was than performed on the detected TEs by considering the three TEs' sub-classes.The distributions of the polarization parameters are shown in Fig. 8.For all the three sub-classes, we observe the prevalence of RL solutions higher than 0.7, thus assuring reliable values of the other polarization parameters.Moreover, the azimuth distributions are in all the cases peaked at about 60 • , which is close to the value found for the class B of series I.In addition, for the sub-class a, another azimuth peak appears at about 150 • .On the other hand, the dip distributions are different for the three sub-classes.Specifically, the sub-class a shows one maximum in the dip distribution at about 10 • -20 • , the sub-class b one maximum at about 20 • and 60 • , and the sub-class c one maximum at about 60 • .Hence, two preferential dip angles are detected, with an increase of the fraction of the lower dips (which correspond to steeper oscillations) moving towards lower frequencies.Main features of TEs in series III are summarized in Table 1.

Modelling TEs
The simplicity of the TEs' waveforms and their frequency content suggest a quite simple source mechanism.Specifically, infrasonic recordings of strombolian-like volcanic explosions show very similar characteristics both in time and frequency domain (e.g.Vergniolle et al., 1996;Johnson et al., 2004).In that case, the source mechanism is associated with the vibrations of gas bubbles before escaping from the vent and bursting.We associate the same source mechanism with TEs.This model can also explain why we observe high-RL polarization, which hence should be associated with P-waves.On the other hand, seismic transient events in "hot" volcanic areas induced by vibrations of fluidfilled cavities such as conduits or dykes (long-period events) display broader-band spectra and longer-lasting waveforms (e.g.Chouet, 2003).Instead, we hypothesize a very efficient acoustic wave/seismic wave coupling, which allows to detect bubble vibrations in the seismic signals and preserve the source-induced polarization.
We will model the TEs as the effect of resonating spherical gas bubbles.Spherical bubbles in an infinite liquid have 3 different modes of oscillation involving shape and volume (Lu et al., 1989) that are responsible for acoustic or seismic emissions.The simplest mode of vibration involves volumetric variations only, leaving the spherical geometry of the bubble unchanged.In this case, restoring force is due to compressibility of gas in the bubble.Source appears as an isotropic point source (unipolar source) and produces a harmonic damped oscillation.Assuming that isothermal heat transfer occurs in the gas at pressure P 0 , that oscillations are small, and long wavelength approximation holds (kR 1 where k is the seismic/acoustic emission wavenumber and R is the bubble radius), the fundamental resonance frequency f 0 of the bubble is given in the approximate form: where γ is the ratio of specific heats of the bubble gas (1.32 for methane), P 0 the ambient pressure (in the assumption that bubble is in equilibrium with hosting fluid) and ρ the ambient density of fluid (of the order of 2000 kg m −3 for mud following Brown, 1990).In the case that pressure equates lithostatic/hydrostatic load, one has where g is the gravitational acceleration and h is the depth where the bubble resonance occurs.In the case of methane bubbles one obtains The most interesting feature of Eq. ( 1) is that resonance frequency inversely depends on the bubble radius: this implies that spectral signatures can be used to evaluate dimension of the resonant bubbles.Excitation of such vibration mode may be obtained in various ways such as, for example, by the detachment of a bubble from a seeping conduit (Manasseh et al., 2001).
Other vibration modes exist (shape modes) that are controlled by surface tension (e.g.Leighton, 1994) mode (Lu et al., 1989).When bubbles interact with fluid/gas interfaces, further vibration modes occur and bubble vibrations become more complex, including bursting and generation of interface gravity waves (Leighton, 1994;Lu et al., 1989;Spiel, 1992).However, when floating bubbles are considered, these radiating modes are much less efficient source of seismic vibrations in the fluid where bubble floats, despite of the fact that their excitation is easier than in the case of volumetric mode (Lu et al., 1989).
The above results can be used to evaluate the dimension of the bubbles responsible for the observed TEs.For this purpose, the depth of the "bubbling reservoir" must be fixed to apply (Eq.3).Following considerations provided by Hovland et al. (1997) and Kopf et al. (2010), this reservoir is considered as the unique source of fluids at the mud cones and salsa emissions.This implies that all the series are recording the effects of a common source.It is in line with the dip solutions from polarization (which point towards deep directions) and with the detection of a common class of TEs for the series I and series III.
A first constraint on the source depth can be obtained considering that no TEs were detected at the salsa (Fig. 1).This could be the effect of attenuation of the seismic signal preventing detectability of the bubbling signal relatively far away from the source.The fact that TEs are only observed at the mud cones implies that their source should be located relatively far away from the salsa and near the volcano summit, possibly immediately below the mud cones area at a depth much shallower than the distance between this seepage structures and salsa (about 300 m).This idea is also in line with the model provided by Roberts et al. (2011) who locate the main vent of the volcano just below the summit area.
In order to fix an upper bound of the source depth, we have at first estimated the mean amplitude of the TEs and of the background tremor.Then, under the hypothesis that TEs are not detected at salsa due to attenuation effects and that they are generated just below the main crater, we have adopted a model for the amplitude decrease taking into account spreading and scattering effects.Specifically, we have adopted A(x) = K 1 (1/x) exp(−K 2 x), where K 2 has been estimated by fixing the wave frequency, the P-wave velocity and the medium quality factor, whereas K 1 was varying on a grid.In this way, we have observed that A (300 m) is equal to the mean amplitude of the background noise when the explosions' mean amplitude is assumed from A for x = 50 m, which we assume to be the maximum source depth.On the other hand, if the reservoir also feeds water in the salsa, it should not be shallower than 10-15 m below the summit area (i.e. the difference between gryphons and salsa heights) to prevent water gushing that is not observed.This lower bound is also confirmed by the polarization analysis that consistently provides relatively steep emerging directions for most TEs.
The fact that gas is ultimately produced in deep seated sediments (few km from the surface in the case of Dashgil vol-cano) does not contrast with the presence of a shallow secondary reservoir where bubbling may occur.In the case of a well-studied mud volcano in Italy (Nirano, Northern Italy), detailed geophysical prospecting revealed the presence of such a reservoir 25-30 m below the main vents (Accaino et al., 2007).
If one assumes a depth in the range 10-50 m, Eq. ( 3) provides bubble sizes of the order of 1-2 m and 0.5-1.5 m respectively for the resonance frequencies of the Class A (3 Hz) and B (5 Hz).These sizes are in line with preliminary estimates (Albarello, 2005) based on theoretical considerations suggesting that a lower bound of the order of 0.1 m exists for the radius of methane carrier bubbles uprising in a mud volcano conduit.An indirect confirmation of these estimates also comes from observed gas bubbles at some mud cones in Dashgil summit area with sizes ranging from few cm up to several tens of cm (Hovland et al., 1997;Etiope et al., 2004b).

Estimates of methane seepage
From data in Table 1, one can estimate an average gas flow from seismic measurements by associating TE frequency and bubble radius via Eq.( 3).Under the hypotheses that (1) the average rate of events detected is representative of all the bubble production at depth; (2) each TE is related to a single different bubble; (3) gas seepage from depth only occurs by bubbling (due to very low permeability of mud); and (4) our experimental observations are representative of the mean behaviour of the volcano, one can estimate the overall mean volume rate produced at depth.As one can see in Table 1, large uncertainty affects bubble radius that is reflected in a larger uncertainty (about one order of magnitude) in the volume associated with each bubble.Bearing in mind these limitations, two estimates of the main gas volume carried out by bubbles in 2006 and 2007 were obtained.In particular, given M the number of classes of TEs detected in the time span T and n i the number of TEs with central frequency f i , the respective mean rate φ corresponding to those TEs is provided by the formula where h is the presumed depth of the resonant bubble.Two possible extreme values of h were considered (10 and 50 m) to evaluate lower and upper bounds for φ.By considering the figures provided in the previous sections for the number of TEs in the different classes and sub-classes (each characterized by a given frequency f 0 ) in the two time series (I and III), the mean rate is estimated in the range [0.9-9] × 10 3 m 3 h −1 for measurements carried out in 2006 (series I) and in the range [0.5-6] × 10 2 m 3 h −1 for measurements in 2007 (series III), i.e. one order of magnitude lower than that observed in 2006.These estimates are largely higher than those provided by Kopf et al. (2010) but are in line (at least when 2007 flow rates are considered) with those provided by Etiope et al., 2004b, and the very preliminary ones by Guliyev and Feizullayev (1997) and Jakubov et al. (1971).
The difference between measurements carried out in 2006 and 2007 concerns the features of the TEs (3 Hz events disappear in 2007), their rate of occurrence (the number of Class B events is one order of magnitude lower in 2007) and the background tremor level (about 30 % lower in 2007 than in 2006).All these features suggest that some major variation in the mud volcano activity has occurred in between.On the other hand, the presence of large variations in gas seepage was also outlined by surface measurements where variations of one order of magnitude are revealed within few hours (Kopf et al., 2010).Of course no sound conclusion can be drawn since there was no correspondence between surface gas seepage measurements and seismic estimates of bubbling (see the time series in Kopf et al., 2010).The presence of large variations in gas flow was also pointed out by Guliyev and Feizullayev (1997).Larger gas outflow in 2006 could also explain the presence of the 3 Hz TEs that are lacking in 2007 measurements.Indeed when a larger seepage from depth occurs, one can expect that the number of bubbles per unit volume increases (see in Table 1 the larger production rates for class-B TEs in 2006 with respect to 2007) increasing the possibility of bubble coalescence.This could result in formation of larger bubbles responsible for lower-frequency TEs.

Discussion and conclusions
Monitoring the activity of the Dashgil mud volcano (Azerbaijan) in its "dormant" phase was carried out through seismic measurement of the volcanic tremor.Measurements were performed by using a very compact, cheap and easy to handle seismometer at three sites in the volcanic crater in 2006 and 2007, revealing the presence of clear energy bursts (transient events -TEs) very distinct from the background noise.The data analysis revealed that these transients are constituted by two classes of nearly monochromatic, linearly polarized seismic events.They have been revealed in the summit area only, where a set of peculiar seepage structures (mud cones and gryphons) exists.Polarization properties of the TEs and their waveforms and spectra suggest that these events can be considered as the surface effect of resonating bubbles in a shallow reservoir just below the volcano summit.
These findings are in line with the common opinion that bubbling plays a significant role in driving mud volcano activity.Indeed, mud volcano tremor exhibits features similar to "strombolian" activity of "hot" volcanoes revealed by acoustic emissions.Some important differences exist anyway.Seismic emissions relative to strombolian activity of hot volcanoes are characterized by remarkable complexity con-cerning both the spectral features of the intense background tremor and the waveforms of the typical energy bursts.This reflects complex dynamics processes occurring in the active structures, involving high energy gas-rock-magma coupling phenomena (Chouet, 1988;Julian, 1994;Fujita et al., 1995;Garcés et al., 1998;Konstantinuou and Schlindwein, 2003;Balmforth et al., 2005;Palo et al., 2009;Matoza et al., 2010;De Lauro et al., 2011), vibrations of magma-filled cavities, such as volcanic conduits, magma chambers or cracks (Crosson and Bame, 1985;Chouet, 2003;De Lauro et al., 2012), etc.None of these phenomena is expected to occur in mud volcanoes where low energy processes only take place (at least during the quiescent phase) as the effect of a slow gas seepage from depth through a mud-filled conduit.Further indications in this direction are the lack of any evident specific spectral signature in the background noise (see, as an example, Figs. 2 and 4) and the fact that it appears to be completely controlled by meteorological and environmental conditions only (see, for example, the strong effect of wind in series III).
This inherent simplicity is also revealed by the examination of the TEs, which show a very uniform and simple waveform.In general, transient "strombolian" events are characterized by low-frequency vibrations accompanied by smaller, high-frequency emission (Vergniolle and Brandeis, 1996).This pattern has been interpreted as the result of bubble vibrations at the surface of magma in the volcanic conduit.The high-frequency signal is lacking in the TEs detected at the Dashgil volcano.This suggests that probably the resonance of the bubble occurs within the fluid and not at its surface.Another possible interpretation is that high-frequency signals are simply too low to be detected at distance due to low-pass frequency filter effect of material damping.Anyway, all these features and the attitude to not consider too complex models when simpler ones provide satisfactory results ("Occam's razor"), supporting the idea that the simple model we adopted for TEs can be considered as fully satisfactory at least at this stage of our studies.
One could wonder why relatively few kinds of bubbles (size) were actually detected instead of a "distribution" of sizes.This could be the effect of relatively high viscosity of mud, which prevents the rising of bubbles having relatively small dimensions (e.g.Albarello, 2005).These will remain trapped in the mud (at rest in the conduit) until their size reaches any critical dimension (for instance, by coalescence of diffusing particles; e.g.Bottiglieri et al., 2005), allowing their uprise.If this occurs at a relatively shallow depth and if the density of the bubbles is not large enough, coalescence is also limited and this also prevents the formation of larger bubbles or slugs.Alternatively, the preferential bubble size can be induced by material constraints in the plumbing system, such as the roof of a reservoir, which promotes the coalescence up to a critical size dictated by the geometrical characteristics of the "trap" (e.g.Chouet, 2003).The fact that bubbles of different size are observed at vents and

D. Albarello et al.: Monitoring methane emission of mud volcanoes
salsas possibly depends on the fractioning occurring when large bubbles reach small conduits near the surface in the presence of water-rich and low viscosity mud.
A basic limitation of the present study (actually shared by most studies concerning gas outflow from inland mud volcanoes) is that, due to inherent logistic difficulties and lack of strong financial support, single station spot measurements were executed only.In this way it was not possible to identify the area where bubble resonance actually occurs.For this purpose, the deployment of a seismic array is necessary and it will represent a basic tool for future researches.A seismic array will improve the source location and will help to remove from the signals potential path and site effects, which at this stage have not been taken into account.In this regard, our results should be considered as a first step in defining a new approach in the study of the gas seepage at mud volcanoes.Moreover, longer-lasting seismic survey are required to establish the stability of the source process.In the present situation, only reasonable hypotheses can be advanced about the location of the reservoir where TEs are generated.Geomorphological considerations and very rough seismic attenuation estimates strongly suggest that such a reservoir should be shallow (10-50 m of depth) and is possibly located at the top of the main vent of the volcano in correspondence with gryphon area in the mud volcano caldera.Some further (though inconclusive) constraints on the source position of the main gas seepage vent can also be inferred by polarization azimuths.Looking at the TEs of 2006, we found slightly different azimuths for the class A and class B of TEs, which, on the other hand, shared the same dip angle equal to about 60 • .Taking into account the bounds of the source depth discussed above and considering P-waves travelling through a homogeneous medium, we get source locations at 17-86 m from the seismometer.The preferential azimuths moreover indicate that the sources of the TEs are located along the directions individuated by the angles 10 • -30 • (or 190 • -210 • ) with respect to the reference direction of instrument (see the map of Fig. 1).Differently from the data of 2006, the polarization analysis of the TEs of 2007 indicates that the subclasses a, b, and c shared the same most frequent azimuth, which is equal to about 60 • or 240 • , but two preferential dip angles at about 10 • and 60 • appeared.This last finding could indicate the same epicentral distance with different source depths.Conversely, if the source depth is common to all the TEs of 2007, the epicentral distance should be different.In detail, following the source and medium hypotheses adopted above, it would be equal to 3-13 m and 14-70 m, respectively, for dips equal to 10 • and 60 • .Moreover, the source epicentres would be located along the direction individuated by the preferential azimuth, i.e. about 60 • .We remark that the sub-class a displays also another maximum in the azimuth distribution at about 170 • , which is mainly associated with low dip angles (< 30 • ), indicating that a minor contribution to the seepage in 2007 comes also from this direction.
In the case of stable epicentral distance and variable source depths, we can image a sort of rising bubble path with two preferential heights for bubble formation.From dip distributions, we inferred that the lowest dip angles are associated with lower frequencies which, in our model, correspond to larger bubbles.It is in line with a higher viscosity of the mud at higher depth (induced by a lower quantity of water) which promotes the coalescence of gas bubbles and, in turn, the formation of larger gas bubbles enhancing the buoyancy effects.Anyway, our model provides differences in sizes between the gas bubbles of the three sub-classes of few tens of centimetres or less, suggesting slight differences in the source depths among the three sub-classes.
The interpretation of transient events in the seismic tremor in terms of bubble resonance suggests a new approach to estimate gas emissions at the mud volcano.Namely, since transient events are nearly monochromatic, a single bubble radius can be associated with each class of TE.By taking into account uncertainty related to the depth where bubble resonance occurs, one can obtain two orders of magnitude of mean gas emissions during the measuring sessions in 2006 and 2007.The result is of the order of 10 4 and 10 3 m 3 h −1 respectively.These values are in line with the very preliminary estimates provided by Jabukov et al. (1971) and more recently confirmed by direct measurements (Etiope et al., 2004).These values are higher of several orders of magnitude than the estimates provided by Kopf et al. (2010).However, these last measurements concern local emissions at one salsa, which is not located on the main vent (Roberts et al., 2011).Furthermore, it appears that, being the reservoir reasonably located below the volcano summit (i.e. in the mud cones area) and relatively far from the salsa, emission at this site could be considered less representative of gas production at depth.
To some extent, production rates obtained here could represent a lower bound due to the possible incompleteness of TE catalogue deduced from field measurements.However, in our opinion such incompleteness could not change the order of magnitude of our estimates.A severe bias instead could concern the hypothesis that each TE is associated with a different bubble.The fact that most TEs share the same dip and azimuth might suggest that one single bubble could provide many TEs.This hypothesis cannot be simply ruled out but it requires some dynamic conditions responsible for repeated bubble perturbations, leaving the bubble in the same position at depth with the same size (the relevant frequency does not change).Alternatively, stable dip and azimuth angles can be also the effect of the interaction of the gas slugs with stable structures of the feeding systems, such as asperities or inhomogeneities in the shallow plumbing system, as occurs in many "hot" volcanoes during long-period seismic events (e.g.James et al., 2006).In that case, the observed seismic signal would be the effect of the redistribution of the gas flux through these structures, with emission of seismic radiation from them.However, the observed waveforms would imply inefficient resonators and very simple and stable geometries of the vibrating structures.On the other hand, the presence of different classes and sub-classes of events that are present or absent in different phases of activity would imply a multiplicity of feeding systems that activate and de-activate depending on the activity rate.This, of course, cannot be excluded "ex ante" but implies a more complex model for the mud volcano structure that is not strictly requested by data (again the "Occam's razor").Future new measurements eventually carried out in array configuration will allow discriminating the most effective model.

Fig. 1 .
Fig. 1.Simplified map of the Dashgil mud volcano crater field (modified from Mazzini et al., 2009) whose geographic position is shown by the arrow in the inset in the upper right corner (north is upward).Circles indicate positions of ambient vibration measurements carried out in 2006 (I, II) and 2007 (III).Bars inside the circles indicate the reference direction for the 3-D seismometer used for ground vibration monitoring.Arrows indicate the direction of the mud flow.Dashed-dotted A-A line indicates the trace of the section in the inset at the bottom of the figure.

Fig. 2 .
Fig. 2. Samples of the time series of seismic tremor (vertical component) registered at the Dashgil mud volcano.(a) Series I measured on June 2006 (in the small inset aa, a magnification of "transient events" (TEs) observed in the series is shown); (b) Series II (June 2006); (c) Series III (July 2007); magnification of a TE for this last series is reported in the inset cc.In this last inset, light gray and black lines respectively report the original and the filtered (i.e.modulated by RL 4 ) signals.

Fig. 3 .
Fig. 3.Stacking of "events" detected in the Series I (on the left) and corresponding spectra (on the right).The upper plots refer to events of Class A and the lower ones to Class B (see text for details).The events were normalized on the respective amplitude before staking.The spectra refer to the mean pattern obtained for each of the two kinds of event.

Fig. 4 .
Fig. 4. Polarization analysis of the Class A (on the left columns) and Class B (in the right column) detected in Series I. From top to bottom: rectilinearity (RL), polarization azimuth, polarization dip and the shape of the typical event.Each dot in the plots in the first three rows represents results of the polarization analysis obtained by considering a sliding window centred at the respective abscissa for a single segment of the global time series.Time windows of 0.5 s and 0.3 s were considered respectively for the Class A and Class B, with an overlap of 75 %.

Fig. 5 .
Fig. 5. Distribution of the peak frequency relative to the events detected in the Series III after polarization filtering (see the text for details).Three maxima are detected in correspondence with the frequencies 4.5 Hz, 5.5 Hz and 6.5 Hz.

Fig. 6 .
Fig. 6.Typical shape of "events" of the sub-classes a (top), b (middle) and c (bottom) detected in the series III.
the distributions of the maxima of the crosscorrelation function estimated by adopting a master TE of sub-classes a, b, c, respectively are reported.In Fig. 7d-e the same graphs computed by choosing a master event from the TEs of classes A and B, respectively, are plotted.

Fig. 7 .
Fig. 7. (a-c) Distributions of the maxima of the cross-correlation function estimated by adopting a master TE of sub-class a, b, c, respectively.(d-e) Distributions of the maxima of the crosscorrelation function estimated adopting a master TE of class A and B, respectively.

Fig. 8 .
Fig. 8. Distribution of rectilinearity (left column), azimuth (middle column) and dip (right column) of the polarization of events (Subclasses a, b, c respectively) detected in Series III.