Deep long period volcanic earthquakes generated by degassing of volatile-rich basaltic magmas

Deep long-period (DLP) earthquakes observed beneath active volcanoes are sometimes considered as precursors to eruptions. Their origin remains, however, unclear. Here, we present a possible DLP generating mechanism related to the rapid growth of gas bubbles in response to the slow decompression of over-saturated magma. For certain values of the gas and bubble content, the elastic deformation of surrounding rocks forced by the expanding bubbly magma can be fast enough to generate seismic waves. We show that amplitudes and frequencies of DLP earthquakes observed beneath the Klyuchevskoy volcano (Kamchatka, Russia) can be predicted by our model when considering pressure changes of ~107 Pa in a volume of ~103–104 m3 and realistic magma compositions. Our results show importance of the deep degassing in the generation of volcanic seismicity and suggest that the DLP swarms beneath active volcanoes might be related to the pulses of volatile-rich basaltic magmas rising from the mantle.

D eep Long Period (DLP) earthquakes occurring in middle to lower crust and uppermost mantle beneath volcanoes [1][2][3][4][5][6][7][8][9] remain enigmatic and in some cases, are believed to have connection with magmatic activity. Similar to volcanic longperiod (LP) seismicity in general 10 , the DLP earthquake has been considered to be generated by rapid pressure variations within magmatic plumbing systems. Alternatively, the effect of thermal stresses within cooling magma bodies has been considered 11 . The cooling magma stalled beneath the crust can also generate DLP earthquakes by so called "second boiling" or repeated pressurization of volatiles exsolved through crystallization, as has been recently suggested for dormant hot-spot Mauna Kea volcano in Hawaii 9 . However, such cooling-related mechanisms are unlikely for DLP events occurring beneath active volcanoes in association with eruptions. Different possible origins of pressure variations resulting in LP seismicity have been considered 12 including the unsteady magma motion, breaking of mechanical "barriers" 9,13 , rapid degassing, etc. In any case, a reasonable model must provide a physical mechanism generating pressure variation dP(t) consistent with observed seismic waves. This implies that the time scale of these variations must by rather short, i.e., comparable with typical frequencies/periods of seismic waves (e.g.,~1 s). Second condition is that the fluid pressure variation should be strong enough and well coupled with the elastic media. This coupling may imply resonances of fluid-filled cracks or cavities 10 that under certain conditions can result in nearly monochromatic and very long duration signals. At the same time, such "strongly resonant" features are not observed for DLP signals that are characterized by rather short durations.
Here we propose that rapid changes of magmatic pressure near the crust-mantle boundary can be caused by nucleation and growth of gas bubbles in response to the slow decompression of over-saturated magma 14 . A volume of magma saturated with H 2 O-CO 2 volatiles is subjected to slow de-pressurization because of its slow upwelling. This magma first reaches the saturation level and then achieves the critical supersaturation after which the gas bubbles nucleate (Fig. 1a) and grow very fast (Fig. 1b). Fast expansion of the bubbly magma deforms the surrounding rocks which respond elastically on the time scale associated with the bubble growth and magma pressure variations. As a result of this elastic rock deformation, seismic waves are radiated (Additional information provided in Methods) and can be recorded by seismographs installed in vicinity of volcanoes.
The pressure variation in the bubbly magma is simulated using the model that accounts for multiple dissolved volatiles (H 2 O-CO 2 ) and diffusive gas transfer from magma into the growing bubbles. It is based on the full solution of advectiondiffusion equation instead of quasi-static approach that was used before (Additional information provided in Methods) 15 . The bubble growth model is adopted to the case of bubble nucleation in basaltic magma 16 .
We compare the results of our modeling with DLP earthquakes observed beneath the Klyuchevskoy volcanic group (KVG) in Kamchatka, Russia. This volcanic group is one of the largest and most active clusters of subduction-zone volcanoes in the World 17 . KVG eruptions and their precursory periods are accompanied by sustained seismovolcanic activity including volcanic earthquakes 7,18,19 and tremors 20 . We particularly focus on a persistent cluster of DLP earthquakes that occur in a small volume located at~30 km depth beneath the Klyuchevskoy volcano 7,19,21 .
The moment magnitudes (Additional information provided in Methods) of these DLP events range between 1.1 and 2.5 with maximum of their distribution at 1.4 ( Supplementary Fig. 1).
Initial data on volatiles in Klyuchevskoy 22 suggested that primary magmas content 2.2-2.9 wt.% of water. Later, a detailed study of melt inclusions in olivines 23 has shown that parental magma has~3.5 wt% H 2 O and 0.35-0.9 wt% CO 2 . Large increase of water content for some melt inclusions (up to 7 wt% H 2 O) was explained by de-compressional crystallization, accumulation volatiles in the melt phase and consequent slow degassing 23,24 .
However, recent experimental data shows that the volatile content of Klyuchevskoy magma is much larger than the one previously directly measured in melt inclusions due to coupled SiO 2 -H 2 O loss 25 , suggesting that primary magma may contain more than 4 wt% of H 2 O. Single H 2 O volatile phase will result in a small saturation depth, but the addition of~0.6 wt% of CO 2 increases volatile solubility dramatically so that magma becomes supersaturated at pressures of 800 MPa (~30 km depth) that alternatively requires~10 wt% of pure H 2 O.
We perform a parametric study to investigate the influence of volatiles content on the dynamics of bubble nucleation and growth. Our results show that the time scale of the bubble growth is mainly controlled by the gas and bubble content in the magma and under certain conditions can be sufficiently fast to generate seismic waves. In particular, we show that amplitudes and frequency content of DLP earthquakes observed beneath the Klyuchevskoy group of volcanoes can be predicted by our model when considering pressure changes of a few tens of MPa in a volume of 10 3 -10 4 m 3 and magmas containing~4 wt% of H 2 0 and~0.6 wt% of CO 2 . Our results provide evidence for the role of the deep degassing in the generation of long-period volcanic seismicity and suggest that the DLP swarms observed beneath active volcanoes

Results
Volatiles content and the depth of degassing unset. Figure 2 shows how the solubility of the CO 2 -H 2 O mixture varies with pressure. Water and carbon dioxide concentrations were parametrized by polynomial functions of pressure and CO 2 content in the bubble at a fixed temperature of 1230 o C estimated from reversed crystallization of Klyuchevskoy melts from atmospheric conditions to 800 MPa pressure (Additional information provided in Methods) and extrapolated to 1000 MPa (Supplementary Table 1). The saturation point is reached at depths of the crustmantle transition (~30 km; P ≈ 825 MPa) for magmas with the volatile content typical for the Klyuchevskoy volcano, implying that degassing may start at such large depths.
Parameters controlling the time scale of bubble growth.  28 . This means that after nucleation the pressure in the gas bubble will be equal to its saturation value and the melt pressure is lower by ΔP. Due to rapid bubble expansion gas pressure decreases extremely fast while melt pressure starts to increase as the volume of the magma increases. Initially gas pressure drop due to bubble expansion dominates pressure increase due to volatile influx into the growing bubble. After reaching minimum value P g starts to increase, concentration gradients in the melt become smoother and volatile flux decreases. At later stages of the growth the difference between melt and gas pressures becomes small and bubble growth is controlled by the diffusion of volatiles. The water diffusion coefficient is 1-2 orders of magnitude larger than the diffusion coefficient for CO 2 . Thus, larger water content of magma for a fixed pressure require smaller amount of dissolved CO 2 and bubble during growth will suck more H 2 O. Adding more water into initial magma results in H 2 O enriched gas and more vigorous bubble and pressure grows (Fig. 3). The effect of water content is enhanced even stronger in predicted seismograms ( Fig. 3d) with H 2 O depleted magmas resulting in very weak signals. We compare amplitudes of synthetic seismograms computed for a magma source volume of 30,000 m 3 (linear dimension of a few tens of meters) with a real seismogram (Fig. 3e) recorded during DLP earthquake with a magnitude M W ≈ 2 at station LGN located nearly above the source region ( Supplementary Fig. 1). Amplitudes and the frequency content (Fig. 3f) are reasonably well predicted with a model based on 4 wt% water in basaltic magma typical for the Klyuchevskoy volcanic group 25 .
We then perform a sensitivity study of several other parameters on the pressure evolution in the growing bubbles and resulting melt ( Supplementary Fig. 2). Critical supersaturation 28 that is required for bubble nucleation does not change melt-pressure recovery time significantly but will affect the amplitude of the source signal. We consider the melt viscosity range 10-10 5 Pa s 29 . If viscosity is smaller than some threshold its influence on resulting pressure is negligible. Only larger melt viscosities typical for more silica reach melts (10 5 Pa s) introduce some delay in pressure recovery. We assume instantaneous bubble nucleation in the whole batch of magma (Additional information provided in Methods). Thus, the size of the cell from which the bubble is growing is controlled by bubble number density (BND). We consider the BND range 30 between 10 11 and 10 15 m −3 . Increase in BND results in smaller cell sizes as S 0~B ND −1/3 . Melt pressure grows faster for smaller S 0 .

Discussion
While the presented comparison of the observed and modelpredicted seismograms ( Fig. 3d-f) is based on significant simplifications of the source (ignoring realistic geometry and possible resonant behavior 10 ) and the propagation effects (ignoring attenuation and wave scattering 31 ), it shows that the amplitudes and the spectral content of the DLP signals observed at the Klyuchevskoy volcanic group can be explained to the order of magnitude by the bubble nucleation and growth in basaltic magmas according to the performed numerical simulation (Fig. 3d-f). Results of the presented modeling show that in the CO 2 -H 2 O rich basaltic magmas the degassing starts at large depths and is vigorous enough to produce strong and rapid pressure variations that can generate seismic radiation with amplitudes and frequency content comparable with those observed by seismographs during DLP earthquakes. Our results suggest that the DLP swarms observed beneath active volcanoes might be related to the intensification of the deep degassing caused by pulses of fresh CO 2   the same time, magma-cooling-related DLP mechanisms can dominate beneath nearly closed or dormant volcanoes. One of the key features of our model is that the depth of occurrence of DLP earthquakes is related to the CO 2 content in magmas. This is especially interesting considering that global volcanic CO 2 fluxes in modern Earth remain poorly known [32][33][34][35][36][37][38] and are often estimated indirectly based on CO 2 /SO 2 or other ratio proxies, with direct CO 2 observations at volcanoes being technically challenging. Our results suggest that studies of the DLP volcanic seismicity provide additional constraints on the magmatic CO 2 content in the deep roots of volcanoes.

Methods
ρ CO2 ¼ ð0:371 þ 0:13 10 À3 TÞ Á P g þ 1194:65 À 0:4665T; ρ H2O ¼ ð0:22 þ 0:13 10 À3 TÞ Á P g þ 892:2 À 0:357T;   Fig. 3 Modeled dynamics of the bubble grows and magma pressure change. Results are shown for the bubble number density of 10 13 m −3 , four different water contents indicated with wt% values in respective plots, and for CO 2 content computed for 828 MPa (red circles in Fig. 2). a Evolution of the bubble radius. b Evolution of magma pressure P m (P g values are shown with gray lines). c Evolution of the CO 2 content in bubbles. d Ground velocities estimated for a source located at a 30 km distance from the receiver (Additional information provided in Methods). e Example of real seismogram (east-west component at station LGN, Supplementary Fig. 1). f Fourier amplitudes computed from synthetic and real (gray line) signals.
Here t is time, r is the radial coordinate, R is the radius of the bubble, v r is the radial velocity, c c and c w are the mass concentrations of CO 2 and H 2 O in the melt, D c and D w are the volatile diffusion coefficients 39 , P g is the pressure of the gas inside the bubble, P m is the melt pressure, σ is the surface tension, µ is the magma viscosity, S is the radius of the cell, G is the shear modulus of the host rock, ρ g is the density of the gas in the bubble that depends on the pressure, temperature T and bubble volatile composition x b CO2 . The densities of pure CO 2 (ρ C02 )and H 2 O (ρ H2O ) are approximated at a limited P-T range using tables produced by NIST Chemistry WebBook (https://webbook.nist.gov/chemistry/).
Equation (2) is subjected to two boundary conditions: concentration gradients are equal to zero at the outer surfaces of the cell mimicking symmetry of the system. At r = R(t) volatiles in magma are in chemical equilibrium with the bubble. Thus, c s ¼ c eq s p; T; x b CO2 À Á . We use D-compress software 40 in order to calculate equilibrium concentrations.
The nucleation time of bubbles t n from a supersaturated melt is related to the bubble number density BND via the nucleation rate I(m −3 s −1 ): t n = BND/I According to classical nucleation theory 41 , I increases extremely fast with oversaturation pressure ΔP: I $ expðÀ1=ΔP 2 Þ. It depends on the temperature and a number of melt properties including surface tension, volume and concentration of water molecules in the melt, as well as distance between them, diffusion coefficient of volatiles at the bubble-melt interface, probability that a nucleus at the top of the barrier will go on to form the new phase, rather than dissolve (Zeldovich factor), and others. With a huge uncertainty of these parameters and difficulties in their experimental constrain, the estimated nucleation rate values vary by orders of magnitudes. For a basaltic melt with an overpressure about 40 MPa a value of I1 0 26 m −3 s −1 has been suggested 42 . With this I-value, nucleation time for BND = 10 13 m −3 (values preferred in our study) is~10 −13 s, which is many orders of magnitude below the typical time scale of the simulated bubble growth and of the observed periods of seismic waves (~1 s). These estimations were obtained assuming that magma degassing is dominated by homogeneous nucleation. In the presence of crystals, their interfaces serve as a preferable location for the heterogeneous nucleation which takes roughly the same time, but produces significantly lower number of bubbles. In the case of heterogeneous nucleation, a pressure perturbation, induced by a limited number of new created bubbles, propagates through a magma-filled cavity providing a trigger for the homogeneous nucleation in the whole volume of magma. Such combination of heterogeneous and homogeneous nucleations is often assumed for many natural systems 28 . The duration of this process is controlled by a propagation time of a pressure pulse across a volume of over-saturated magma. With typical dimensions of a few tens of meters and sound speed being of the order of a few km/s, the combined heterogeneous and homogeneous nucleation will take less than 0.01 s, i.e., two orders of magnitude below typical bubble growth times. Therefore, we consider instantaneous nucleation in the whole volume.
Numerical method. Equation (1) can be integrated analytically and gives the following velocity distribution in the melt phase: v r ¼ dR dt R 2 r 2 . In order to solve Eq. (2) in a fixed domain we use front-fixing method 43 ð Þ gives extra advective term in Eq. (2). The resulting equation is discretized on an irregular 1D mesh with a decrease of the step size towards the growing bubble boundary (ξ = 0). The resulting system of equations with threediagonal matrix is solved by means of Thomas algorithm 44 . The forward step starts at the outer domain boundary (ξ = 1). The linear relation of volatile concentration on the bubble boundary and in the nearest mesh points together with discretized Eqs. (3) and (4) allows to calculate all parameters on the bubble-melt interface. Then, concentration distribution in the whole domain is calculated during backward substitution. We found this method stable and computationally efficient in comparison with explicit methods that require extremely small timesteps for stability reasons.
Estimation of magma composition. In order to estimate magma compositions in the deep magma reservoir we used "Petrolog" software 45 . Reverse crystallization from a more evolved magma (sample 12KY-108-1, 1987 AD eruption 46 ) was performed. The starting pressure is set to atmospheric level and magma is H 2 Osaturated. Our simulations reveal total amount of mineral phase of 20% for the starting composition, which is in a good agreement with the measurements on the samples 47 . Incremental increase in pressure to 800 MPa leads to the change in composition presented in Supplementary Table 1. These values were obtained considering the volatile component composed only of H 2 O, resulting in a 800 MPa magma containing almost 11 wt% of dissolved water. Adding even a small amount of CO 2 affects significantly the water solubility that can be reduced to a few wt% as shown in Fig. 2 along the 800 MPa isobar. Based on data about Klyuchevskoy magma volatile content [22][23][24][25]48 , we retain for our modeling a composition with~4 wt% of H 2 0 and~0.6 wt% of CO 2 .
Estimation of magnitudes of deep low-frequency earthquakes. The DLP signals whose energy is concentrated in a narrow spectral band between 1 and 2 Hz are dominated by S-waves (Fig. 3e, f). The seismic moment can be approximately estimated from maximal signal amplitude in the following way. We start with an expression of the far-field (hypocenter distances exceeding 10 wavelengths) S-wave displacement 49 u S and ignore the radiation pattern assuming that it approximately averages to 1. Based on this we can relate the time derivative of seismic moment with the observed S-wave displacement: where t is time M 0 is seismic moment, ρ is density, β is S-wave speed, and r is the hypocentral distance. The observed ground velocity v S is the derivative of the displacement that for a nearly monochromatic signal can be approximately estimated via multiplication by 2πf max : where f max is the dominant signal frequency. Integration of Eq. (6) to obtain the whole seismic moment can be also approximately estimated with dividing by 2πf max . This leads to a final expression used to approximately estimate the seismic moment from one station: where v s max is the maximum amplitude of velocity seismograms (taking into account all three components). The final estimate is averaged from several stations that recorded the earthquake. We use f max = 1.5 Hz and typical crustal values for density 27 , ρ = 2830 kg m −3 , and seismic velocity 50,51 , β = 3500 m s −1 . The moment magnitude M W is then computed as: Estimation of seismic radiation emitted by expanding magma volume. For simplicity, we start with considering a volume with a perfectly spherical shape embedded in an infinite elastic space with bulk modulus K. In response to the magma pressure change dP(t), the volume will be modified by dV(t): For a perfectly spherical magma body, the volume change can be related to the seismic moment as 49 : A spherically symmetric source would radiate in the far field only P waves. At the same time, signals from real DLP earthquakes are dominated by S waves. A simple explanation of this observations can be related to the deviation of the magma body shape from a perfect sphere. In this case, the change of the magma pressure will induce a significant amount of shear stress in the surrounding rocks resulting in a strong S-wave radiation 52 . A possible example is a pure tensile crack mechanism for which the seismic moment tensor can be written as 49  where λ and μ are Lamé constants that for most of elastic solids are nearly equal and have the same order of magnitude as bulk modulus (K = λ + 2/3μ) implying that to the order of magnitude the relationship (11) between seismic moment (observed amplitudes of waves), pressure variations, and volume of affected fluid remain valid. Seismic radiation from such source for many directions is dominated by S-waves 53 . At this stage, we do not consider detailed description of seismic radiation from a non-spherical source that would vary significantly depending on the exact magma volume shape. We rather make an order of magnitude estimation and consider that Eq. (11) describes the relationship between the magma pressure change and the seismic moment observed in the far field (hypocenter distances exceeding 10 wavelengths). Based on Eq. (6), the ground displacement can be expressed as: and the ground velocity is computed as its time derivative.

Data availability
The seismological time series used for the analysis were provided by the Kamchatka Branch of the Geophysical Survey of Russian Academy of Sciences (GS RAS) and are available on request (http://www.emsd.ru). The data are not publicly available due to the internal regulation of the GS RAS.