Mathematical model for volcanic harmonic tremors

Harmonic tremors consist in the release of infrasonic energy associated with volcanic activity. The typical frequency range of harmonic tremors is 0.1–12 Hz. We suppose that the harmonic tremors are due to the formation of bubbles entrapped in cavities that oscillate converting thermal energy into mechanic energy. Reproducing the natural phenomenon through an experimental apparatus, we propose here a mathematical model to describe the oscillatory mechanism and to detect the frequency as a function of the main physical parameters. We show that the frequency obtained through the model is in agreement with the one obtained through experimental measurements and with the data available from the literature, proving the consistency of the proposed model.

Harmonic tremors with frequency 0.1-12 Hz often occur in volcanos related to eruptive phenomena, but they can also be present in volcanoes without any unusual activity 1 . Infrasounds are commonly reported as time series or spectrograms, to show the evolving frequency content in a signal 2,3 . Some of these tremors are attributed to being the result of more and more rapidly occurring pulses, while others are the result of an explosion. Finally, some harmonic tremors just emerge and disappear apparently not having a known cause.
The measurement of acoustic (infrasound) signals is important to localize the source of tremors and to gather a deeper insight on the phenomena responsible of their emergence 4 . The source is localized using a grid searching analysis as reported in 4 . To understand the dynamics of volcanic tremors, the fumarole thermal trends can also be studied. Thermal measurements are generally done by using infrared thermometers, and by investigating the non-linear behaviour in the temperatures we can discover the link between multiple fumaroles and the effect of the fluid dynamics at depth.
Recently Busse et al. 5 have proposed that harmonic tremors are linked to thermoacoustic instabilities, i.e. the oscillations of large gas bubbles entrapped in cracks. This association is quite straightforward because volcanoes are surrounded by hydrothermal reservoirs, namely rock cracks in which vapor or gas bubble are heated by the magma. We however have to remark that the causes of the harmonic tremors are still uncertain and that a better understanding of possible sources would be a further step in expanding our knowledge [6][7][8][9][10] .
There are many non-linear behaviors possible in a volcanic system related to both fluid flow and chemical system, as commented e.g. in 11 . During the fast uprising of the gas, steam, liquid droplets and solid particles the gas temperature may span in the interval 80-1100 °C and turbulent and laminar motion alternate, generating oscillations related to fluid dynamic features and posing a challenge to the fluid sampling [11][12][13][14][15] .
A recent investigation 16 carried out at Pisciarelli fumarole, Campi Flegrei volcano (Italy), reports the data from a seismic station located at 8 m distance from the fumarole. Polarization analysis results indicate that the signal is dominated by the vertical component, and originates from the vent. Its prominent spectral peak at ~10 Hz frequency is stable over time. A characteristic RSAM (Real-time Seismic-Amplitude Measurement) value is computed daily by assuming as representative the minimum value registered during the night. The daily characteristic RSAM values (in m s −1 ) are used by the authors to obtain the reduced average seismogram of the tremor time series. This reduced signal (known as RSAM or fumarolic tremor) includes the tremors generated by the entire fumarolic system at Pisciarelli, the subterraneous feeding channels and the surface vents. During the period of observation, the amplitude of fumarolic tremor at Pisciarelli increased by more than one order of magnitude with many peaks that, according to the authors, are linked to the hydrothermal fumarolic activity and reflects the increasing temperature/pressure conditions of the hydrothermal system.
Other experimental studies carried out at Stromboli (Italy) 17,18 , where continuous (strombolian) volcanic activity is present, measured gas velocity at fumaroles outlet of 10 m/s and pressure and thermal oscillations with frequency of ~2 Hz and average amplitude of 20 Pa and 30 °C.
In this paper we propose that a possible source of volcanic tremors is the oscillations of vapor bubbles confined within cavities of the hydrothermal basin surrounding the volcano. The idea is not new (see 1 ) and can explain phenomena occurring at shallow depths (i.e. not exceeding few km from the ground surface), so that steam can separate from the liquid. The rocks hosting an hydrothermal system show a large amount of fractures of millimeters and sub-millimeters width and length from few centimeters up to meters (e.g. 19 ) in which vapor bubbles can form. Bubbles spontaneous oscillations induced by heat flux have been a research topic in the last two centuries (see 20 ). In the hydrothermal systems the heat source could be the hot fluid upraising in a fumarolic duct (e.g. Pisciarelli), or the magma that seeps in cavity, fractures and volcanic conduits. The magma, while equilibrating its temperature with the one of the host rocks, releases the heat needed to support the bubbles oscillations, that efficiently convert the thermal energy into mechanical energy. Since volcanic tremors of frequency 2-10 Hz are a very common phenomenon observed worldwide, they can be related to fundamental physical effects such as bubble oscillations. The general idea is that vapor bubbles develop near the heat source and, under proper conditions, dissipate energy not only by conduction but also by performing mechanical work (i.e. oscillating).
The dynamics of a single vapor bubble in an incompressible or nearly incompressible liquid is a famous problem. In the authors knowledge, the first theoretical model was introduced in 21 , who considered an incompressible fluid. The majority of these works is reviewed in 22 . We also acknowledge the study of 23 devoted to the collapse of a spherical vapor bubble in a subcooled liquid, and the one of 24 in wich the growth and collapse of a vapor bubble inside a microtube is experimentally and theoretically studied. Vapor bubbles dynamics in microchannels is considered also in 25 . The problem considers a heater element producing vapour and a cooler at the interface. The solution of this free boundary problem is found for a steady-state condition in a rectangular geometry. A number of studies have dealt with the growth or collapse of an axisymmetric bubble near a plane wall 26 and with the instability problem of an oscillating spherical bubble 22 . Many important experimental and theoretical studies have investigated the problem of motion and deformation of cavitation bubbles in water [27][28][29][30][31] . The results obtained so far, show that the asymmetry produced for large bubbles (due to the shear stress close to the wall) affects the dynamic evolution of the bubble itself.
In this study we mimic the natural phenomenon using an experimental apparatus and then we develop a mathematical model capable of describing the observations. The experimental set-up reproduces in the lab the bubble oscillations occurring in a rock cavity. The experiment, described in Section 1, consists in a capillary tube filled with pure water and heated from the bottom. The lateral wall of tube is properly insulated, so that the heat flow occurs only from the bottom of the tube. When the water reaches the boiling temperature we observe the formation of a vapor bubble which grows up to a stable configuration around which it performs oscillations of small amplitude. We focus on this stage recording, by means of infrasonic detection system, the sound produced by the bubble oscillations. The analysis of the frequency spectrum revealed a peak in the spectrogram around 12 Hz.
To provide an insight on the bubble dynamics, a mathematical model of the system is presented in section 2. The aim of the model is to give a simple but rigorous description of the oscillatory behavior of this simple system that acts as a proxy for bubbles in fractures under a confining pressure (represented by the water column). The model is based on the energy balance between the liquid and vapor phase, where the vapor is treated as an ideal gas. To describe the phase transition between the two phases of the same single constituent (water) we have used the Clausius-Clapeyron relation. The FFT (Fast Fourier Transform) of the oscillating signal provided by the model exhibits a peak around 12 Hz, showing an excellent agreement with the spectrum obtained from the experimental measurements.

the experiment
The simple experiment proposed here, carried out in the IGG laboratories, evidences the dynamic of a vapor bubble on a lab scale. We designed this experiment to be as simple as possible using: (i) An amorphous silica tube 1 m long, 8 mm diameter, wall thickness 1 mm. The amorphous silica tube is needed for safety, because it is able to resist to both the thermal stress and the shock induced by bubble vibrations; (ii) Pure water (MilliQ grade); (iii) Pyro-Bloc modules (Lynn Manufacturing inc.); (iv) Bunsen burner (air-methane).
(v) Nikon D3000 buit-in condenser microphone. (vi) OriginLab software for the determination of the spectrogram of the recorded signal.
The diameter of the silica tube is not very important as long as the tube diameter is small enough to prevent convective motion but sufficiently large to minimize the meniscus effect. During our experiments, we tested different pipes (6, 8, 10, 12 mm of diameter) and in the 12 mm tube diameter the phenomenon do not present, likely due the bubbles upraising and generating convective cells.
The tube is suspended over a Bunsen burner and then inserted in a 5 cm thick PyroBloc module where a hole, having the same diameter of the tube, has been previously drilled. The tube is heated from the bottom since the 5 cm thickness PyroBloc practically avoids any lateral heating. The silica tube is filled with pure water (MilliQ grade) up to nearly 60 cm. After having set up the system, the Bunsen burner is turned on and positioned with the oxidizing flame at the tube bottom (Figs 1 and 2).
After approximately 30 minutes a vapor bubble of nearly 3 cm height is formed at the bottom of the tube. For about 30 minutes after the bubble forming, the system performs small oscillations. The vapor bubble acts www.nature.com/scientificreports www.nature.com/scientificreports/ as a "spring" for the liquid column above it, which performs oscillations with amplitude of the order of a few millimeters. Then the bubble collapses and, in a relatively short time, a second bubble which grows very rapidly is formed. The expansion of the second bubble eventually produces a sort of "explosion" that makes the water column outpour from the tube.
The experimental data shown here refer to the second stage, that we named "oscillatory stage". In Fig. 3 the spectrogram of the recorded periodic signal is reported, showing a peak around 12 Hz.

the Mathematical Model
We consider a one dimensional model as the one depicted in Fig. 4, and focus on the second stage (i.e. the "oscillatory stage"), so that the initial time t = 0, corresponds to the time at which the vapor bubble has just formed. We denote by L the height of the liquid column and by y the vertical coordinate (pointing upward) so that y = 0 is the tube bottom and y = h(t) represents the height of the vapor bubble. We assume that in the vapor phase, [0, h], the physical variables do not depend on y and model the liquid column, [h, h + L], as a body moving with uniform velocity. The mathematical problem is formulated writing the energy balance for the liquid and vapor phase (coupled with appropriate interface conditions) and the motion equation for liquid column. We remark that on the interface h(t) (which is a free boundary) a vapor-liquid phase transitions takes place, so that h(t) is not a material interface.
where θ a is the ambient temperature, θ(t) is the temperature of the vapor in the bubble and θ in (y) is the temperature profile at t = 0. We remark that (2) 2 provides the coupling between the temperature field in the water and the the bubble vapor temperature whose evolution equation will be illustrated in Section 2.2. Introducing the transformation and the new variable Equation (1) and conditions (2) become www.nature.com/scientificreports www.nature.com/scientificreports/ the vapor bubble. We model the vapor as an ideal gas in which temperature and pressure are linked by Clayperon's equation 32 . Hence we write where p is the vapor pressure, ρ is the vapor density, λ is the latent heat of vaporization, r is the specific gas constant and θ o is the evaporation temperature at the reference pressure p o . From (6) we see that given by where δQ is the infinitesimal increment of heat supplied to the system and pdV is the infinitesimal work done by the system, V being the volume of the bubble. We notice that V(t) = Ah(t), where A is the cross section area of the capillary tube. In our case, since the lateral wall is adiabatic and no heat source/sink is present, we have where  Q b > 0 is the heat supplied at the bottom surface by the bunsen burner (constant in time) and  Q h is the heat flux through the interface y = h(t). The infinitesimal variation of internal energy dU in the vapor phase is given by dU = d(ρVc v θ), where c v is the specific heat at constant volume. Recalling that c p − c v = r (c p being the s p ecific heat at constant pressure) (9) rewrites as Exploiting (6) and (8) we write so that (10) becomes where ρV is the mass of the vapor within the bubble. To close Eq. (11) we need an expression for  Q h . To this aim, we recall that h(t) is an evaporation/condensation interface: here the heat fluxes have a jump proportional to the heat absorbed, or released, by the phase change. Applying the Stefan condition (see 33 where  Q l is the heat flux of the liquid phase at the interface y = h which, exploiting Fourier's law, is given by with θ l obtained solving (3), (4). In conclusion Eq. (11) can be rewritten in the following form l l is the mass of the liquid column, p a is the atmospheric pressure acting on the free surface of the liquid column and where the last term represents the damping due to viscous friction. The right hand side of (14) is therefore the net force acting on the water column. Choosing the reference pressure as (the reference pressure is thus the sum of the ambient pressure and the pressure exerted by the water column) with H ≈ 10 m, and recalling (6), we may rewrite (14) as the non dimensional problem. Let us rescale the problem with non dimensional variables. We introduce where t c is a characteristic time still to be selected and h o is the characteristic height of the bubble. Next we set with p o reference pressure given by (16), θ o the water boiling temperature at pressure p o and  where t D = L 2 D −1 is the characteristic diffusive time in the water. The boundary and initial conditions become We then we introduce where t b is the characteristic time of the heat supply due to the bunsen burner and t h is the characteristic time of phase transition, so that (13) becomes so that (21) can be rewritten as To write the non dimensional version of the momentum balance of the water column we introduce the characteristic oscillation time  Since we are interested in the oscillation of the bubble we select t o as characteristic time, that is t c = t o . We end up with the following mathematical problem for the unknowns (We omit the "~" to keep notation light) Θ(ξ, t), ζ(t) and h(t) where the values of Table 2 are determined using the physical parameters of Table 1. The initial and boundary conditions are www.nature.com/scientificreports www.nature.com/scientificreports/ System (25) is quite involved because two nonlinear Ordinary Differential Equations (ODEs) are coupled with a parabolic equation for Θ. We however remark that ν − 10 9 , and so only a thin layer of water facing the interface is heated by the vapor. This allows to approximate the heat flux by (see Appendix A) where ξ d is the rescaled thickness of the thermal boundary layer in which temperature decreases from the boiling temperature to ambient temperature. Setting v =  h, the system (25) can be simplified getting to the following first-order system of ODE's where ρ(ζ) is given by (22) and F(ζ, t) is given by (27).

numerical Results
Here we plot the solution of system (28) using the parameters of Tables 1 and 2 and using the following initial data In particular we plot the height h(t) and the temperature θ(t) of the bubble as a function of time, see Figs 5 and 6. Finally we perform a FFT (Fast Fourier Transform) of the oscillation signals and evaluate the frequency domain. The power spectrum is shown in Fig. 7 where we can see a peak at ~12 Hz, i.e. at the same frequency recorded during the experiment (see Fig. 3). In particular, Fig. 6 shows that the bubble temperature is essentially constant during the oscillations, meaning that the phase change process has little influence on the dynamics. So, we can consider o o with ρ o given by (18) and o o p o being the "equilibrium" pressure given by (16). This fact allows to simplify considerably the mechanical model (14). Indeed, neglecting the damping term, Eqs (14), (6) 2 , (29) and (30)   www.nature.com/scientificreports www.nature.com/scientificreports/ Assuming that  y 1 we approximate the above with whose frequency of oscillation is Referring to the data of Table 1 we obtain ν ≈ 12 Hz, i.e. the recorded frequency and the output of the complete model (28). In particular, (32) leads to some of observations that we shall discuss in the next section.

Discussion and Conclusions
The formation of bubbles in a liquid is an efficient way of mechanical energy generation and this fact has been recognized in geophysics 35 . Boiling of groundwater causes the formation and growth of bubbles in the liquid permeating the rocks of the geothermal and volcanic system. The bubbles start to oscillate if supported by a suitable heat flux. In particular, the bubbles oscillation occurs at a frequency which is derived neglecting both the phase change process and the temperature variations. So, from (14) we derive (32). Equation (31) can be related to Mathieu equation 36 that is known to produce resonance. In our interpretation, among the many possible www.nature.com/scientificreports www.nature.com/scientificreports/ oscillating frequencies occurring in the natural system, it is likely that a number of different bubbles begin to resonate producing a harmonic tremor at a given frequency.
As one can easily recognize from (32), the frequency is a function of the characteristic bubble height and the column length. A 3D plot of the function (32) is shown in Fig. 8. The pair (h o , L) spans in the domain [0.01, 1] × [0. 5,20] meaning that we are considering 1 cm as the minimum height of the bubble and 1 m as the maximum, while minimum L is 0.5 m and maximum L is 20 m. As one can see, the frequency range is [0.6, 12.09] which is exactly the range of the frequencies for harmonic tremors recorded all over the world. To provide a better understanding on how we may accurately model the measured frequency for bubble oscillation and the effect of H/L in Eq. (32), we model the fumarolic tremor measured by 16 at Pisciarelli fumarole, where the frequency at fumarole head is nearly 10 Hz. We model the system as our lab experiments, using the same conditions (the real conditions are very similar, only the temperature is higher) with a H/L ratio of 0.056667. Setting Q h = 0.1 W (see Table 1), we are able to reproduce the measured frequency span.
The heat flux affects the timespan of the phenomenon, meaning that Q h = 10 −3 W (see Table 1) is sufficient to sustain the oscillation, while with larger values of Q h the bubble size increases, modifying the oscillation frequency, until it disappears. This effect provides a frequency span for the oscillation of a single bubble, that in the natural system could overlap with the resonance between different bubbles oscillating at different frequencies. At Vulcano, Vulcano Island (Italy) 37 many seismic signals related to bubbles in the hydrothermal system are detected. In particular, these signals are produced by fluid flow inside cracks and conduits with frequencies that depend on the gas volume fraction. Among them there are monochromatic events that show a single frequency peak at 6 or 8 Hz and generally last about 10-15 s, with few of them lasting up to 30-50 s with a slow amplitude decay. According to the author, these features are highly similar to the so-called "tornillos" found in other volcanoes (Galeras 38 , Tongariro 39 , Tatun Volcano Group 40 ).
To model the 6-8 Hz oscillations we need a H/L (parameters from (32) and Table 1) ratio of 0.05, with a heat flux Q h of 0.1 W. If we use a H/L of 0.0125 with an heat flux Q h of 0.1 (more than 50 s) or 0.4 to reduce the lifetime of the phenomenon, we obtain oscillations in a frequency range of 2-4 Hz corresponding to the higher infrasonic signal frequency detected in 37 . Ripepe et al. 4 measured infrasonic waves with similar spectral content between 2 and 10 Hz at Stromboli volcano crater. The peaks centered at 5 Hz can be attributed to gas bubbles expansion, that is the phenomenon studied in this paper. Therefore we may state that our model is capable of reproducing the measured infrasound signals for the case study of Solfatara (Pisciarelli fumarole), Vulcano (Vulcano Island) and Stromboli. In particular, the case of Vulcano covers typical infrasound signals found in many other volcanoes around the world. In this work, we are able to reproduce the measured infrasound signals from different volcanoes (namely Solfatara at Pisciarelli fumarole, Vulcano, Stromboli) in which the infrasound signals are attributed to bubble oscillations in cracks or in magma (Stromboli). We remind that one of the main challenge is to prove the relation between the emitted infrasound and the bubble oscillations, but once this is done we are able to accurately reproduce the signal frequency with our model.
where we recall that Θ a is the rescaled ambient temperature. The above means that in a layer of thickness ξ d the temperature of the liquid is the boiling temperature, while outside the layer the fluid is everywhere at ambient temperature. A reasonable value for ξ d is 1/60. Indeed, recalling that the length of the column is 60 cm, the assumption ξ d = 1/60 means that we have 1 cm of the water column close to the bubble at the boiling temperature, while the rest is essentially at ambient temperature. Introducing the new non dimensional variables , a we get (see Table 2)