Effect of heterogeneities on evaluating earthquake triggering of volcanic eruptions

Recent researches have indicated coupling between volcanic eruptions and earthquakes. Some of them calculated static stress transfer in subsurface induced by the occurrences of earthquakes. Most of their analyses ignored the spatial heterogeneity in subsurface, or only took into account the rigidity layering in the crust. On the other hand, a smaller scale heterogeneity of around hundreds of meters has been suggested by geophysical investigations. It is difficult to reflect that kind of heterogeneity in analysis models because accurate distributions of fluctuation are not well understood in many cases. Thus, the effect of the ignorance of the smaller scale heterogeneity on evaluating the earthquake triggering of volcanic eruptions is also not well understood. In the present study, we investigate the influence of the assumption of homogeneity on evaluating earthquake triggering of volcanic eruptions using finite element simulations. The crust is treated as a stochastic media with different heterogeneous parameters (correlation length and magnitude of velocity perturbation) in our simulations. We adopt exponential and von Karman functions as spatial auto-correlation functions (ACF). In all our simulation results, the ignorance of the smaller scale heterogeneity leads to underestimation of the failure pressure around a chamber wall, which relates to dyke initiation. The magnitude of the velocity perturbation has a larger effect on the tensile failure at the chamber wall than the difference of the ACF and the correlation length. The maximum effect on the failure pressure in all our simulations is about twice larger than that in the homogeneous case. This indicates that the estimation of the earthquake triggering due to static stress transfer should take account of the heterogeneity of around hundreds of meters.


Introduction
A causal relationship between large earthquakes and following volcanic eruptions has been indicated by many researches (e.g.Linde and Sacks, 1998;Nostro et al., 1998;Hill et al., 2002;Manga and Brodsky, 2006;Battaglia et al., 2012).It is important to understand the triggering mechanisms quantitatively not only for the scientific interests but also the mitigation of volcanic risks.Since a volcanic eruption requires a change in stress for permitting the migration of magma to the surface, the stress change induced by an earthquake is considered to be a key issue of the triggering mechanism.Diez et al. (2005) calculated the static stress change following three Mw 5.2 earthquakes occurred near Cerro Negro volcano, Nicaragua, and indicated that the reduction in minimum horizontal principal stress was sufficient to trigger the eruption.Walter and Amelung (2007) calculated volumetric strain induced by historical large earthquakes at subduction zones, i.e. the 1952 Kamchatka M 9.0, the 1960 Chile M 9.5, the 1964 Alaska M 9.2, and the 2004 Sumatra-Andaman M 9.3.They found that the eruption rate increased in the region of expansion associated with the megathrust earthquakes.Wang et al. (2011) calculated the changes of stressstrain field resulted from the 2011 Tohoku-oki earthquake M 9.0 using a 3-D linear elastic half-space (Okada, 1992), and pointed out that Mt.Fuji, Shinmoedake and Changbaishan volcanoes, which started volcanic unrest or eruptions immediately after the earthquake, are located in regions of volumetric expansion.Feuillet et al. (2011) investigated the effect of earthquakes on the volcanoes in the Lesser Antilles considering Coulomb stress, normal stress and pressure changes.Their results showed that the most recent intraplate earthquakes may have brought the increment of the static stress Published by Copernicus Publications on behalf of the European Geosciences Union.
beneath the volcanoes pumping system and subsequent eruptions and seismic activities.
In their analyses, the numerical models were treated as homogeneous media.However, seismic surveys and geological investigations have indicated that the subsurface include various scales of heterogeneities which might have the influence on the estimation of the stress fields in the subsurface generated by earthquakes.Zhao et al. (2004) investigated deformation and stress generated by strike-slip, tensile and thrust faulting, considering the rigidity layering in the crust.They showed that the surface deformation and stress changes are distorted by the ignorance of crustal structural information.Furthermore, an oversimplification of the hypothesis of the homogeneity in volcanic areas, which might lead to misinterpretations, has been indicated in not only earthquake triggering studies (e.g.Manconi et al., 2010).
On the other hand, the smaller spatial scale heterogeneity of around hundreds of meters in the crust has been indicated (e.g.Line et al., 1998;Chung et al., 2009).Different from the larger scale heterogeneity like the rigidity layering, it is difficult to introduce that kind of heterogeneity in analysis models accurately because velocity fluctuations tend to be random at a scale lower than hundreds of meters.In other words, even if we can understand the heterogeneity parameters (e.g. the correlation length) in the target area, it is difficult to specify the concrete distributions of them, especially in the deeper crust shaded by the complex scattering surface structure.Most of analyses have neglected that kind of heterogeneity in their numerical models.However, the smaller scale heterogeneity might also perturb the stress field induced by earthquakes around volcanic area similar to the disturbance from the larger scale heterogeneity like the rigidity layering.Thus the stress disturbance due to the smaller scale heterogeneity could influence the earthquake triggering.
In the present study, we investigate the effect of the ignorance of the smaller spatial scale heterogeneity on evaluating the earthquake triggering of volcanic eruptions by means of calculating failure pressure change around a magma chamber wall (Pinel and Jaupart, 2003;Albino et al., 2010), using finite element simulations.We conduct numerical simulations with different heterogeneous parameters, i.e. the difference of the auto-correlation functions (ACF), the correlation length and the degree of the heterogeneity, and investigate whether the ignorance of the smaller spatial scale heterogeneity is an oversimplification or not.

Model descriptions
Our two-dimensional numerical model is shown in Fig. 1.Host rock is treated as a stochastic heterogeneous medium.The mechanical properties of the host rock are characterized by P-wave velocity V p , mass density ρ and Poisson's ratio σ .P-wave velocity is represented as the superposition of a background component V p0 and a random stochastic component δV p (x).The stochastic component is represented by its spatial ACF, which is commonly Gaussian, exponential or von Karman (Frankel and Clayton, 1986).We adopt exponential and von Karman functions, whose ACFs are where the correlation length L and the magnitude of the velocity perturbation ε. ε is defined as ε 2 = δV p (x) 2 /V 2 p0 .ν is Hurst number, ( ) is Gamma function and K ν ( ) is modified Bessel function of order ν.Mass density and Poisson's ratio are fixed to 2700 kg m −3 and 0.25, respectively.Circular magma chamber, whose diameter D is fixed to 800 m, is embedded in the host rock.We also introduce a strong heterogeneous region near the chamber with diameter D .The bulk modulus K of magma depends strongly on content exsolved gas (e.g.Tait et al., 1989).In this study, the bulk modulus of magma is fixed to 10 GPa.
The governing equation in this study is a static equilibrium equation as follows: where ∇ • is a divergence operator, s is the stress tensor and f is the external force.The stress tensor is represented as follows: where C and γ are the elastic and strain tensors, respectively.Elements of C can be calculated by V p , ρ and σ .Equation ( 3) is discretized by a finite element method (FEM).The normal displacement at the lower boundary of the model is fixed to zero.Enforced deformation is applied to the lateral edges as a large scale strain caused by an earthquake.In the present study, we consider volumetric expansion which is induced by megathrust earthquakes in the area where volcanoes erupted (Walter and Amelung, 2007; Wang et al., 2011).The large scale strain is fixed to 1.0 × 10 −6 .Although using a different strain value affects the magnitude of the calculated stress field, the relative magnitude is not changed.We consider brittle failure at the magma chamber wall for the initiation of dykes which causes magma migration towards the surface.Following Pinel andJaupart (2003, 2005) and Albino et al. (2010), we calculate the deviatoric component of the minimum compressive stress at the chamber wall Pr and pressure change in the magma chamber Pc.Comparing these two pressure changes, we can investigate the magma chamber moves away from failure conditions or not.
Positive and negative values of the difference between Pr and Pc indicate that the inhibition and approach of dyke initiation, respectively (Albino et al., 2010).
We change the correlation length, the magnitude of the velocity perturbation, ACFs and the diameter of the strong heterogeneous region in order to investigate the effect of the smaller scale heterogeneities on evaluating the earthquake triggering of volcanic eruptions.

Correlation length and magnitude of velocity perturbation
In this section, we conduct finite element simulations without the strong heterogeneous region near the chamber, i.e.D = 0 m, in order to investigate the effect of the correlation length and the magnitude of the velocity perturbation.We use the exponential function (Eq. 1) as the ACF in this section.The effect of the choice of the ACF is investigated in the next section.Figure 2 shows examples of the P-wave velocity distributions of our model with different L and ε.
In order to investigate the effect of the correlation length, the ratio R between L and D (R = L/D) is varied by only changing L.
Figure 3 shows the results of our simulations for the failure pressure change (Pr − Pc) with different velocity perturbations.The ratio R is fixed to 1 (i.e.L = 800 m).The error bar denotes the standard deviation for 20 model realizations.We can observe that the heterogeneous host rocks always induce the smaller failure pressure than that from the homogeneous host rock.This indicates that the assumption of the homogeneity could underestimate the failure pressure at the chamber wall.Not only the average of the failure pressure but also the standard deviation of that increases with increasing the magnitude of the velocity perturbation.This means that the larger magnitude of the velocity perturbation possesses great  Figure 4 shows the results with different correlation lengths.The magnitude of the velocity perturbation is fixed to 10 %.There is little change in the magnitude of the failure pressure although the standard deviation slightly increases with increasing the ratio R.This means that the larger spatial scale of the heterogeneity has a possibility to induce the lower failure pressure.Nonetheless, the correlation length has a smaller effect on the failure pressure at the chamber wall than the degree of the heterogeneity.
We conducted further simulations under different combinations of the heterogeneous parameters, L and ε, however, the above trend does not change in any case.

ACF
We next investigate the effect of the ACFs on the failure pressure.In this section, we adopt a von Karman function (Eq.2) as the ACF.We change the Hurst number ν in Eq. ( 2) from 0 to 0.5.The von Karman functions with the Hurst number of 0 and 0.5 correspond to the self-similar and exponential functions, respectively.The magnitude of the velocity perturbation ε, the correlation length L and the diameter D are fixed to 10 %, 800 and 0 m, respectively.Figure 5 shows examples of the velocity distributions with different Hurst numbers from the same random number seed.The model with the smaller Hurst number includes the shorter-scale fluctuations.Figure 6 shows the results with the different Hurst numbers.In each case, the averages and the standard deviations are not so different between each case.As shown in Fig. 4, the correlation length has an insignificant effect on the failure pressure.The results of Fig. 6 agree with that of Fig. 4, i.e. the same magnitude of the velocity perturbation brings the similar failure pressure.

Change of the heterogeneity near the chamber
In this section, we investigate the effect of a strong change of the heterogeneity near the chamber.We change the diameter of the strongly heterogeneous region from 0.9 to 1.5 km.The correlation length is fixed to 800 m.The magnitude of the velocity perturbation near the chamber is set to twice larger than the back ground region, which is fixed to 5 %.Thus, the magnitude of the velocity perturbation near the chamber is 10 %.The exponential function is used as the ACF.An example of the velocity distribution with D = 1.5 km is shown in Fig. 7.
Figure 8 shows the results with the different diameters of the strongly heterogeneous region near the chamber.The horizontal axis is the normalized diameter (i.e.D /D).As shown in the figure, the diameter D has an insignificant effect on the failure pressure.The averages and the standard deviations are similar to the case with ε = 10 % in Fig. 3.This means that the heterogeneity in the immediate vicinity of the chamber is a dominant factor rather than the heterogeneity of the background host rock.

Discussions and conclusions
We investigated the effect of the assumption of the homogeneity in the subsurface on the failure pressure at the magma chamber wall due to the occurrence of earthquakes using an FEM.Our results showed that the ignorance of the heterogeneity, whose spatial scale is around hundreds of meters, has a certain degree of effect on the estimation of the tensile rupture at the chamber wall induced by a megathrust earthquake.It is important to note that the ignorance of the heterogeneity always brings underestimation of the triggering effect regardless of the heterogeneous parameters.The stress field perturbation due to the heterogeneity makes a local stress concentration around the chamber wall, whereas the stress relaxation area is also brought.Since the failure occurs at the closest part to the critical state, the heterogeneity could play an important role on evaluating the earthquake triggering of volcanic eruptions.
We investigated the effect of the correlation length, the degree of the heterogeneity, the ACFs and the strong change of the heterogeneity near the chamber on the failure pressure at the magma chamber wall induced by an earthquake.Among them, the degree of the heterogeneity has the largest effect, i.e. the larger magnitude of the velocity perturbation leads to the smaller failure pressure.On the other hand, the correlation length and the difference of the ACF have almost no effect on the failure pressure.The diameter of the strong heterogeneous region near the chamber also has an insignificant effect, whereas the magnitude of the heterogeneity around the chamber strongly affects the failure pressure.This indicates that the magnitude of the velocity heterogeneity only around the magma chamber is a dominant factor.
The maximum effect of the heterogeneity in the failure pressure for all our simulations is about twice larger than the homogeneous case when the magnitude of the velocity perturbation ε is 20 %.De Siena et al. ( 2011) investigated the scattering attenuation structure below Campi Flegrei Caldera using the autocorrelation functions.In their analysis, the space averaged value of ε 2 was estimated to 0.06 ± 0.02.This value is higher than the magnitude of the heterogeneity used in our numerical simulations.This indicates that the greater effect on the estimated failure pressure may be possible and the heterogeneity should be taken into account in stress transfer analyses of the earthquake triggering of volcanic eruptions.

Fig. 1 .
Fig. 1.The schematic figure of two-dimensional finite element simulations.A circular magma chamber is embedded into the model at a depth of 3.2 km.A shaded region is a strong heterogeneous region near the chamber.The model is elongated by the extensional displacement as a large scale strain induced by a megathrust earthquake."Rollers" at the bottom boundary indicate zero normal displacements.

Fig. 3 .
Fig. 3.The failure pressure with the different velocity perturbations ε.The vertical and horizontal axes represent the failure pressure and the velocity perturbation, respectively.The unit of the vertical axis is Pa.Filled circles and error bars are average and standard deviation for 20 model realizations.Dotted horizontal line represents the failure pressure in the homogeneous case.

Fig. 4 .
Fig. 4. The failure pressure with different correlation length.The horizontal axis represents the logarithmic ratio between correlation length and diameter of the magma chamber.The other descriptions are same as Fig. 3.

Fig. 6 .
Fig. 6.The failure pressure with different Hurst numbers.The horizontal axis represents the Hurst number ν in Eq. (2).The other descriptions are same as Fig. 3.

Fig. 7 .
Fig. 7. Example of the P-wave velocity distribution with D = 1.5 km.The unit of the contour is m s −1 .L and ε are 800 m and 10 %, respectively.

Fig. 8 .
Fig. 8.The failure pressure with different diameter D .The horizontal axis represents normalized diameter D /D.The other descriptions are same as Fig. 3.