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.

Introduction

Deep Long Period (DLP) earthquakes occurring in middle to lower crust and uppermost mantle beneath volcanoes[1-9] remain enigmatic and in some cases, are believed to have connection with magmatic activity. Similar to volcanic long-period (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 H2O–CO2 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.

Fig. 1

Conceptual model of fluid-related source of long-period earthquakes.

a Bubble nucleation in a volume of magma saturated with H2O–CO2 volatiles. b Bubble and pressure growth deforming the surrounding rocks.

The pressure variation in the bubbly magma is simulated using the model that accounts for multiple dissolved volatiles (H2O–CO2) and diffusive gas transfer from magma into the growing bubbles. It is based on the full solution of advection-diffusion 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% H2O and 0.35–0.9 wt% CO2. Large increase of water content for some melt inclusions (up to 7 wt% H2O) 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 SiO2-H2O loss[25], suggesting that primary magma may contain more than 4 wt% of H2O. Single H2O volatile phase will result in a small saturation depth, but the addition of ~0.6 wt% of CO2 increases volatile solubility dramatically so that magma becomes supersaturated at pressures of 800 MPa (~30 km depth) that alternatively requires ~10 wt% of pure H2O.

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 ~103–104 m3 and magmas containing ~4 wt% of H20 and ~0.6 wt% of CO2. 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 might be related to the pulses of fresh CO2–H2O rich basaltic magmas rising from the mantle.

Results

Volatiles content and the depth of degassing unset

Figure 2 shows how the solubility of the CO2–H2O mixture varies with pressure. Water and carbon dioxide concentrations were parametrized by polynomial functions of pressure and CO2 content in the bubble at a fixed temperature of 1230 oC 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 crust-mantle 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.

Fig. 2

Gas saturation isobars as function of CO2–H2O content.

Thin blue lines show saturation isobars for different pressures (indicated values in MPa). Thick solid line indicates the decompression path of the Klyuchevskoy magmas[24] from initial state at 1 GPa shown with a star. Red circles show compositions along the 828 MPa isobar with 2, 3, 4, and 5 wt% of H2O tested with numerical modeling (results shown in Fig. 3).

Fig. 3

Modeled dynamics of the bubble grows and magma pressure change.

Results are shown for the bubble number density of 1013 m−3, four different water contents indicated with wt% values in respective plots, and for CO2 content computed for 828 MPa (red circles in Fig. 2). a Evolution of the bubble radius. b Evolution of magma pressure P (P values are shown with gray lines). c Evolution of the CO2 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.

Parameters controlling the time scale of bubble growth

Figure 3a, b shows typical evolution of bubble size, gas and melt pressure in the basaltic magma with density[26] of 2800 kg m−3, viscosity of 10 Pa s and containing 1013 bubbles m−3 and for four different concentrations of H2O in initial magma. The CO2 contents (red circles in Fig. 2) were computed for initial pressure of 828 MPa that corresponds to the lithostatic pressure at depth of ~30 km for the average crustal density[27] of 2830 kg m−3. Based on the experimental observations we adopt that the critical supersaturation for the bubble nucleation corresponds to the over-pressure of ΔP = 40 MPa[14,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 Pg 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 CO2. Thus, larger water content of magma for a fixed pressure require smaller amount of dissolved CO2 and bubble during growth will suck more H2O. Adding more water into initial magma results in H2O 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 H2O depleted magmas resulting in very weak signals. We compare amplitudes of synthetic seismograms computed for a magma source volume of 30,000 m3 (linear dimension of a few tens of meters) with a real seismogram (Fig. 3e) recorded during DLP earthquake with a magnitude MW ≈ 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–105 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 (105 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 1011 and 1015 m−3. Increase in BND results in smaller cell sizes as S0 ~ BND−1/3. Melt pressure grows faster for smaller S0.

Discussion

While the presented comparison of the observed and model-predicted 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 CO2–H2O 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 CO2–H2O rich basaltic magmas rising from the mantle. This mechanism supports that the DLP earthquakes are early seismic manifestations of activation of deep parts of the Klyuchevskoy volcano plumbing systems. Similar behavior might be expected in other open and very active volcanic systems (with adjusting the model parameters based on their magma compositions and volatile contents). At 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 CO2 content in magmas. This is especially interesting considering that global volcanic CO2 fluxes in modern Earth remain poorly known[32-38] and are often estimated indirectly based on CO2/SO2 or other ratio proxies, with direct CO2 observations at volcanoes being technically challenging. Our results suggest that studies of the DLP volcanic seismicity provide additional constraints on the magmatic CO2 content in the deep roots of volcanoes.

Methods

Mathematical model of gas bubble growth

We consider growth of an individual bubble in the center of a spherical cell of melt that expands with the bubble and supports it with volatiles. The spherically symmetric model includes equations of mass conservation of the melt in a cell Eq. (1), diffusion equations for volatiles (H2O–CO2) Eq. (2), Rayleigh–Lamb equation for bubble growth with negligibly small inertia terms and the equation for the melt pressure evolution due to expansion of the surrounding elastic host rock Eq. (3), mass balances for volatiles in the bubble Eq. (4), and equations that describe physical properties of the components Eq. (5):Here t is time, r is the radial coordinate, R is the radius of the bubble, v is the radial velocity, c and c are the mass concentrations of CO2 and H2O in the melt, D and D are the volatile diffusion coefficients[39], P is the pressure of the gas inside the bubble, P 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, ρ is the density of the gas in the bubble that depends on the pressure, temperature T and bubble volatile composition . The densities of pure CO2 (ρC02)and H2O (ρ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, . We use D-compress software[40] in order to calculate equilibrium concentrations.

The nucleation time of bubbles t from a supersaturated melt is related to the bubble number density BND via the nucleation rate I(m−3 s−1): t = BND/I According to classical nucleation theory[41], I increases extremely fast with oversaturation pressure ΔP:. 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 I~ 1026 m−3s−1 has been suggested[42]. With this I-value, nucleation time for BND = 1013 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: . In order to solve Eq. (2) in a fixed domain we use front-fixing method[43]. A coordinate transformation 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 three-diagonal 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 H2O-saturated. 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 H2O, resulting in a 800 MPa magma containing almost 11 wt% of dissolved water. Adding even a small amount of CO2 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-25,48], we retain for our modeling a composition with ~4 wt% of H20 and ~0.6 wt% of CO2.

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]uS 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 M0 is seismic moment, ρ is density, β is S-wave speed, and r is the hypocentral distance. The observed ground velocity vS is the derivative of the displacement that for a nearly monochromatic signal can be approximately estimated via multiplication by 2πfmax:where fmax is the dominant signal frequency. Integration of Eq. (6) to obtain the whole seismic moment can be also approximately estimated with dividing by 2πfmax. This leads to a final expression used to approximately estimate the seismic moment from one station:where 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 fmax = 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 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.