Finite element analysis of ground deformation due to dike intrusion with applications to Mt. Etna volcano

A 2D finite elements study was carried out to analyse the effects caused by dike intrusion inside a heterogeneous medium and with a realistic topography of Mt. Etna volcano. Firstly, the method (dimension domain, elements type) was calibrated using plane strain models in elastic half-spaces; the results were compared with those obtained from analytical dislocation models. Then the effects caused both by the topographic variations and the presence of multi-layered medium on the surface, were studied. In particular, an application was then considered to Mt. Etna by taking into account the real topography and the stratification deduced from seismic tomography. In these conditions, the effects expected by the dike, employed to model the 2001 eruption under simple elastic half-space medium conditions, were computed, showing that topography is extremely important, at least in the near field. Mailing address: Dr. Rosario Occhipinti Amato, Dipartimento di Ingegneria Industriale e Meccanica, Facoltà di Ingegneria, Università degli Studi di Catania, Viale A. Doria 6, 95125 Catania, Italy; e-mail: occhsaro@hotmail.com


Introduction
Volcanic areas are zones where marked ground deformation often occurs in response to the action of volcanic sources and magma uprising mechanisms.A classic mechanism of magma uprising is related to intrusive processes characterising the emplacement of dikes.A dike can be represented by an opening crack.Tensile dislocation theory has been applied with success to model the case of a rectangular shaped source (fig.1a) with tensile opening in an elastic homogeneous half-space (Okada, 1985;Yang and Davis, 1986).The application of dislocation theory to the case of tensile mechanisms has often been used in the last ten years and applied successfully to model ground deformations by inverting the recorded data to infer source parameters (e.g., Yang et al., 1992;Bonaccorso, 1999).However, the real characteristics of the medium are not perfectly homogeneous.Furthermore, real topography is very different from the flat reference surface assumed in analytical solutions (fig.1b), whereas an analytical solution in realistic situations with topography and stratifications of the medium is impraticable.A valid alternative is represented by the use of the Finite Elements Method (FEM).The study was carried out to investigate the effect of a tensile crack in a volcanic edifice, through FEM numeric technique.The use of this analysis has been consolidated in other fields such as engineering, physics, etc., while its application is more recent in geophysics and volcanology (e.g., Dieterich and Decker, 1975;Melosh and Raefsky, 1981).In Finite element analysis of ground deformation due to dike intrusion with applications to Mt. Etna volcano general, boundary conditions, optimal size of the domain, meshing, characteristics of the medium, and the type of the elements to be used are not known a priori, so it is necessary to calibrate the models.Thus, a robust preliminary set-up is needed to obtain realistic results.A finite elements analysis is always preceded by the meshing of the domain.This is a very delicate phase that influences the accuracy of the calculation, since an erroneous choice of the domain size and the type of elements would lead to unreliable results.The present study initially focussed on the analysis of a very simple two-dimensional domain, with the aim of testing the choice of elements.Computations are carried out with the FE software MARC (MARC, 1994).The results were compared with the theoretical ones expected from the dislocation theory in a homogeneous medium (e.g., Okada, 1985).Then two-dimensional analyses were done by introducing in the model the topography of Mt.Etna, along an eastwest section, to examine the effect of the topography on the deformation field.In addition to the topography, the effects of the material and stratification of the medium were studied considering an elastic multi-layered medium with the Young modulus deduced from recent seismic tomography.
The aim was to compare the pattern of deformation obtained analytically from simple half-space medium with that produced by the same source considered inside a realistic layered medium with topography.

Model without topography in a homogeneous medium
As a first step, we considered two 2D FEM models reproducing a rectangular tensile dislocation in a homogeneous half-space.The geometry of the intrusion is shown in fig.1a and in table I.In agreement with other authors (e.g., Chevallier and Verwoerd, 1988), zero-stress was assumed at the upper surface, while zero displacement values were assumed at bottom and lateral boundaries.Material behaviour was assumed elastic with λ = µ, namely Poisson ratio ν = 0.25, and Young modulus of 40 GPa (Paul et al., 1987;Russo et al., 1997).After a number of simulations, a domain of 1 0 0 × 21 km 2 , meshed with triangular elements, proved sufficient to represent the problem under examination (table I).Results highlighted that the  vertical and horizontal deformation fields practically coincide with the analytical ones (e.g., Okada, 1985) (fig.2).This allowed us to verify the validity of the size domain and element type used in this simulation.

Model with topography in a homogeneous medium
We performed numerical modelling with real topography of an East-West section of the Etna volcano (fig.3).This topography permits results from the FEM analysis and the analytical model to be compared (Bonaccorso et al., 2002) obtained by inverting tilt and GPS data recorded during the emplacement and final uprising of the dike of the Etna 2001 eruption.The final parameters of this model are reported in table I.The dislocation model assumed the reference plane at 1500 m a.s.l., namely the mean elevation of points belonging to the geodetic permanent networks.
After a number of calibrations of the models, a domain of 20 ×60 km 2 was considered suffi- cient (table II).Here too, the medium was assumed homogeneous with λ=µ and Young's module equal to 40 GPa.As in the case without topography, the value of Young's module, in the condition of λ=µ, does not alter the values of the surface deformations.Instead, the results show an evident influence of the topography on the deformation fields (fig.4).In particular, the Fig. 6.Comparison between the deformation field of the tensile crack in a homogeneous medium (i.e.without stratifications) with topography and one in a stratified medium with topography.The top of the crack is positioned at 1420 m a.s.l., the geometric parameters are in table I.

1545
Finite element analysis of ground deformation due to dike intrusion with applications to Mt. Etna volcano modulus, we have applied the following equation E=5/6ρV 2 p where Vp is the seismic P-wave propagation velocity in an elastic medium, with ν = 0.25.Results compared to non-stratified medium models (Section 2.2) revealed no significant differences (fig.6).This implies that, as a first estimation in an intrusion model with rectangular source and constant opening, the effects produced in a homogeneous half-space loss of symmetry of the distribution of the vertical and horizontal displacements is recognised, owing to the non-symmetry of the domain.With respect to the half-space model, one can observe a re-distribution of the deformation with higher mean values having a shift of the maximum peaks.In particular, this displacement, also present for the vertical changes, causes a reverse tilt near the zone above the crack, in agreement with other authors (Trasatti et al., 2003).Furthermore, the maximum amplitude of the displacements undergoes an increase in the part where the topography has a wide depression in the surface profile.

Models with topography in a layered medium
Starting from the models with topography, another analysis was performed.We added the stratification of the medium by introducing seven layers of elastic linear materials with different Young moduli (fig.5).The thickness and the elastic moduli (table II) were carried out from a recent seismic tomography of Mt.Etna (Chiarabba et al., 2000).In order to obtain the Young with λ = µ and in a half-space with parallel stratifications with λi = µi, are the same.

Conclusions
Finite elements analysis of a geophysical problem, such as an intrusive mechanism, is complex and therefore requires various simplifications.In 2D models the necessary parameters for the successive simulations were tested: material, surrounding conditions and typology of elements.The results showed a good correspondence with the values obtained from the dislocation theory, which allows us to affirm that the finite elements method provides reliable results if the model is suitably calibrated.If we consider the tensile crack with a constant opening, in the condition λ = µ, the effects of the value of Young's module are practically negligible on the behaviour of the surface deformations.We have shown that the effects of a simple model of tensile dislocation with a uniform constant opening in an infinite single-layer with λ=µ coincide with those caused by the same source in a multi-layered stratification with λi= µi.Instead, the effect of the topography is important, above all in the near-field of accentuated topographical variations, usually present in summit areas.The deformations in the far-field, for instance on the volcano slopes, where the topography has less than 20°slopes, provide comparable values to those expected also in more realistic situations with stratifications and topography.Therefore, the much faster parametric inversions associated with analytical models, which assume half-space conditions, provide reliable preliminary parameters if near-field data are excluded from the inversion.

Fig. 3 .
Fig. 3. Map of Mt.Etna with near E-W section topography and dike position.

Fig. 4 .
Fig. 4. Comparison between the deformation fields of a model in a homogeneous half-space medium without topography and a model homogeneous half-space with real topography.The top of the crack is located at 1420 m a.s.l.; the zero of the abscissa corresponds to the crack projection onto the surface.

Fig. 5 .
Fig. 5. Stratifications of the elastic medium and grid used in the FE calculations; zoom of the summit zone is reported.

Table II .
Parameters of FEM analysis of a dike intrusion for the different medium analysed: elastic homogeneous half-space without topography; elastic homogeneous half-space with topography; medium with topography and seven stratification layers.