Annals 47, 5, 2004 1541 ANNALS OF GEOPHYSICS, VOL. 47, N. 5, October 2004 Key words FEM – ground deformation – dike in- trusion – real medium 1. Introduction Volcanic areas are zones where marked ground deformation often occurs in response to the action of volcanic sources and magma up- rising mechanisms. A classic mechanism of magma uprising is related to intrusive process- es 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 charac- teristics of the medium are not perfectly homo- geneous. Furthermore, real topography is very different from the flat reference surface as- sumed in analytical solutions (fig. 1b), whereas an analytical solution in realistic situations with topography and stratifications of the medi- um is impraticable. A valid alternative is repre- sented by the use of the Finite Elements Method (FEM). The study was carried out to investigate the effect of a tensile crack in a vol- canic 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 geo- physics 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 Rosario Occhipinti Amato (1), Marco Elia (1), Alessandro Bonaccorso (2) and Guido La Rosa (1) (1) Dipartimento di Ingegneria Industriale e Meccanica, Facoltà di Ingegneria, Università degli Studi di Catania, Italy (2) Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Catania, Italy Abstract 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 con- ditions, 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, Dipar- timento di Ingegneria Industriale e Meccanica, Facoltà di Ingegneria, Università degli Studi di Catania, Viale A. Do- ria 6, 95125 Catania, Italy; e-mail: occhsaro@hotmail.com 1542 Rosario Occhipinti Amato, Marco Elia, Alessandro Bonaccorso and Guido La Rosa 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 prelimi- nary 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 ini- tially focussed on the analysis of a very simple two-dimensional domain, with the aim of test- ing 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 dis- location theory in a homogeneous medium (e.g., Okada, 1985). Then two-dimensional analyses were done by introducing in the mod- el the topography of Mt. Etna, along an east- west section, to examine the effect of the to- pography 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 de- formation obtained analytically from simple half-space medium with that produced by the same source considered inside a realistic lay- ered medium with topography. 2. Materials and method 2.1. Model without topography in a homogeneous medium As a first step, we considered two 2D FEM models reproducing a rectangular tensile dislo- cation in a homogeneous half-space. The geom- etry 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 ra- tio ν = 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 km2, meshed with triangular elements, proved sufficient to represent the problem under exam- ination (table I). Results highlighted that the Fig. 1a,b. a) Analytic representation: rectangular tensile dislocation in an elastic homogeneous half-space; b) representation with Finite Elements Method (FEM): intrusion in a homogeneous medium with topography. Table I. Parameters of the analysed dike. This geom- etry is the one of the analytical models obtained by Bonaccorso et al. (2002) to model ground deforma- tion changes recorded by permanent geodetic net- works during the dike emplacement of the Etna 2001 eruption. Dislocation Width: Inclination Top crack opening: B → W angle: δ depth: H [m] [m] [degree] [m] 3.5 2300 90° 80 a b 1543 Finite element analysis of ground deformation due to dike intrusion with applications to Mt. Etna volcano Fig. 2. Comparison between the theoretical deformation field and that found with the FEM model for a tensile dislocation in a homogeneous and elastic half-space. Dots indicate the FE solution, solid lines the analytical solu- tion (grey colour for horizontal displacement and black colour for vertical displacement). The agreement between the two sets of solution is evident. For the geometrical parameters see table I. Fig. 3. Map of Mt. Etna with near E-W section to- pography and dike position. vertical and horizontal deformation fields prac- tically 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. 2.2. 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 analyti- cal model to be compared (Bonaccorso et al., 2002) obtained by inverting tilt and GPS data recorded during the emplacement and final up- rising 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 geo- detic permanent networks. After a number of calibrations of the models, a domain of 20 × 60 km2 was considered suffi- 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. cient (table II). Here too, the medium was as- sumed homogeneous with λ = µ and Young’s module equal to 40 GPa. As in the case without topography, the value of Young’s module, in the 1544 Rosario Occhipinti Amato, Marco Elia, Alessandro Bonaccorso and Guido La Rosa Table II. Parameters of FEM analysis of a dike intrusion for the different medium analysed: elastic homoge- neous half-space without topography; elastic homogeneous half-space with topography; medium with topogra- phy and seven stratification layers. Model Width × No. No. Mean nodes Element Crack Dislocation No. Young Depth × Depth elements nodes distance in type W × H opening: layers module layer the apex area [km] [m] [m] [m] [Mpa] [km] Homogeneous 21 ×100 26 850 13804 4 triangular 2300 ×80 3.5 1 40 000 21 000 without topography Homogeneous 20 × 60 48 788 24 775 20 triangular 2300 ×80 3.5 1 40 000 60 000 with topography Stratification 20 × 60 48 788 24 775 20 triangular 2300 × 80 3.5 7 6800 2000 with 8350 1000 topography 10100 0 12 400 −1000 15 300 −2000 18 700 −3000 22 100 >−4000 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 posi- tioned 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 equa- tion E=5/6ρV 2p 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 sig- nificant differences (fig. 6). This implies that, as a first estimation in an intrusion model with rec- tangular source and constant opening, the ef- fects produced in a homogeneous half-space Fig. 5. Stratifications of the elastic medium and grid used in the FE calculations; zoom of the summit zone is reported. loss of symmetry of the distribution of the verti- cal 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). Fur- thermore, the maximum amplitude of the dis- placements undergoes an increase in the part where the topography has a wide depression in the surface profile. 2.3. 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 sev- en layers of elastic linear materials with differ- ent Young moduli (fig. 5). The thickness and the elastic moduli (table II) were carried out from a recent seismic tomography of Mt. Etna (Chia- rabba et al., 2000). In order to obtain the Young 1546 Rosario Occhipinti Amato, Marco Elia, Alessandro Bonaccorso and Guido La Rosa with λ = µ and in a half-space with parallel strat- ifications with λi = µi, are the same. 3. Conclusions Finite elements analysis of a geophysical problem, such as an intrusive mechanism, is complex and therefore requires various simpli- fications. In 2D models the necessary parame- ters for the successive simulations were tested: material, surrounding conditions and typology of elements. The results showed a good corre- spondence with the values obtained from the dislocation theory, which allows us to affirm that the finite elements method provides reli- able 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 de- formations. We have shown that the effects of a simple model of tensile dislocation with a uni- form constant opening in an infinite single-lay- er with λ = µ coincide with those caused by the same source in a multi-layered stratification with λi = µi. Instead, the effect of the topogra- phy 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 stratifica- tions and topography. Therefore, the much faster parametric inversions associated with analytical models, which assume half-space conditions, provide reliable preliminary pa- rameters if near-field data are excluded from the inversion. REFERENCES BONACCORSO, A. (1999): The March 1981 Mount Etna in- ferred through ground deformation modelling, Phys. Earth Plan. Inter., 112, 125-136. BONACCORSO, A., M. ALOISI and M. MATTIA (2002): Dike emplacement forerunning the Etna July 2001 eruption modelled through continuous tilt and GPS data, J. Geophys. Res., 29 (13), doi: 10.1029/2001GL014397. CHEVALLIER, L. and W.J. VERWOERD (1988): A numerical model for the mechanical behaviour of intraplate vol- canoes, J. Geophys. Res., 93 (B5), 4182-4198. CHIARABBA, C., A. AMATO, F. BARBERI and E. BOSCHI (2000): Recent seismicity and tomographic modeling of the Mount Etna plumbing system, J. Geophys. Res., 105, 10923-10938. DIETERICH, J.H. and R.W. DECKER (1975): Finite Element Modeling of surface deformation associated with vol- canism, J. Geophys. Res., 80, 4094-4102. MELOSH, H.J. and A. RAEFSKY (1981): Simple and efficient method for introducing faults into finite element com- putations, Bull. Seismol. Soc. Am., 71, 1391-1400. MARC (1994): Analysis Research Corporation. User’s Guides, Version 2000 (MARC Analysis Research Corporation, Palo Alto, CA), A-F vols. OKADA, Y. (1985): Surface deformations due to shear and tensile faults in a half-space, Bull. Seismol. Soc. Am., 75, 1135-1154. PAUL, A., J.P. GRATIER and J. BOURDON (1997): Numerical model for simulating deformation of Mount St. Helens volcano, J. Geophys. Res., 92 (B10), 10,299-10,312. RUSSO, G., G. GILBERTI and G. SARTORIS (1997): Numeri- cal modeling of surface deformation and mechanical stability of Vesuvius volcano, Italy, J. Geophys. Res., 102 (B11), 24,785-24,800. TRASATTI, E., C. GIUNCHI and M. BONAFEDE (2003): Ef- fects of topography and rheological layering on ground deformation in volcanic regions, J. Volcanol. Geotherm. Res., 122 (1/2), 89-110. YANG, X. and P.M. DAVIS (1986): Deformation due to a rectangular tension crack in an elastic half-space, Bull. Seismol. Soc. Am., 76, 865-881. YANG, X., P.M. DAVIS, P.T. DELANEY and A.T. OKUMARA (1992): Geodetic analysis of dike intrusion and mo- tion of the magma reservoir beneath the summit of Kilauea volcano, Hawaii, 1970-1985, J. Geophys. Res., 97, 3305-3324. (received July 4, 2003; accepted December 19, 2003)