untitled 119 ANNALS OF GEOPHYSICS, VOL. 51, N. 1, February 2008 Key words Campi Flegrei – deformation – gravity – finite element method – heterogeneous medium 1. Introduction Volcanic activity produces deformation and gravity changes that may be used as precursors of future eruptions. Monitoring of active areas, together with interpretation of these changes, may reveal the physics and characteristics of the magma reservoir. Indeed, the detection of ground deformation will constrain the location and type of deep reservoirs, while gravity mon- itoring is recognized as a valuable tool for map- ping subsurface density distributions. Both of them contribute to quantify the change in sub- surface mass. Several models have been devel- oped to interpret geodetic and gravity signals in active volcanic areas. Walsh and Rice (1979) developed a method to calculate the gravity changes as due to subsurface mass redistribu- tion with respect to the gravimeters at surface. The calculation is suitable for any dislocation source in the half-space; however the stress change in a dislocation-free half-space must be known. Bonafede and Mazzanti (1998) sug- gested the separation of density variations af- fecting the computation of gravity changes into 3 different contributions: the first contribution comes from the input of new mass from remote distance, the remaining two from the linearized continuity equation of the material already present in the region. Fernandez and Rundle Gravity changes due to overpressure sources in 3D heterogeneous media: application to Campi Flegrei caldera, Italy Elisa Trasatti (1) and Maurizio Bonafede (2) (1) Istituto Nazionale di Geofisica e Vulcanologia, Roma, Italy (2) Università degli Studi di Bologna, Bologna, Italy Abstract Employing a 3D finite element method, we develop an algorithm to calculate gravity changes due to pressurized sources of any shape in elastic and inelastic heterogeneous media. We consider different source models, such as sphere, spheroid and sill, dilating in elastic media (homogeneous and heterogeneous) and in elasto-plastic me- dia. The models are oriented to reproduce the gravity changes and the surface deformation observed at Campi Flegrei caldera (Italy), during the 1982-1984 unrest episode. The source shape and the characteristics of the medium have great influence on the calculated gravity changes, leading to very different values for the source densities. Indeed, the gravity residual strongly depends upon the shape of the source. Non negligible contribu- tions also come from density and rigidity heterogeneities within the medium. Furthermore, if the caldera is elas- to-plastic, the resulting gravity changes exhibit a pattern similar to that provided by a low effective rigidity. Even if the variation of the source volumes is quite similar for most of the models considered, the density inferred for the source ranges from ∼400 kg/m3 (super critical water) to ∼3300 kg/m3 (higher than trachytic basalts), with drastically different implications for risk assessment. Mailing address: Dr. Elisa Trasatti,Istituto Nazionale di Geofisica e Vulcanologia, Via di Vigna Murata 605, 00143 Roma, Italy; e-mail: elisa.trasatti@ingv.it Vol51,1,2008_DelNegro 5-02-2009 10:28 Pagina 119 120 E. Trasatti and M. Bonafede (1994) developed a method to compute gravity variations due to point magma intrusion in a horizontally layered elastic-gravitational me- dia. The solutions were extended and general- ized to point-sources in multi-layered vis- coelastic-gravitational media (Fernandez et al., 2001). A comprehensive tool for evaluating de- formation and gravity changes due to disloca- tion in horizontally layered viscoelastic media was also provided recently by Wang et al. (2006). The approaches described above are on- ly partially able to simulate the complex geolo- gy of volcanic regions, since they model point- sources in horizontally layered media. Further- more, biases in estimate of source density may arise from oversimplified models. Trasatti et al. (2005) developed finite element deformation models taking into account the structural and rheological characteristics of Campi Flegrei. The introduction of rigidity layering and plastic rheology in the local crust causes a shrinking of the deformation shape and enhances the maxi- mum values within the caldera. This behavior allows us to estimate a source depth of 5 km (plastic medium) instead of 3 km (elastic medi- um), adopting the same overpressure and in ac- cordance with information from groundwater studies and petrological data. Recently, Curren- ti et al. (2007) calculated the gravity changes in FE models characterized by an ellipsoidal source expanding in a heterogeneous medium with real topography of Mt. Etna. They found that gravity estimates may be biased in terms of mass gain/loss if medium complexities are neg- lected. In the present paper we attempt to model both the deformation and the gravity. We devel- op a technique to calculate gravity variations due to deep pressurized sources in heteroge- neous media. The method consists in two sepa- rated steps: i) the displacement and strain fields due to the source are computed by FE model- Fig. 1. Map of Campi Flegrei caldera showing the outer rim (open triangles) and the inner caldera (full trian- gles), the leveling lines (dotted) were surveyed during 1982-1984 unrest episode, while the squares represents the gravity stations active in that period. Vol51,1,2008_DelNegro 5-02-2009 10:28 Pagina 120 121 Gravity changes due to overpressure sources in 3D heterogeneous media: application to Campi Flegrei caldera, Italy ling; ii) the results are employed in a procedure to integrate numerically the gravity variations from the displacements and strain fields. The great advantage of this technique is that the FE models may include several complexities: arbi- trary source shape, lateral heterogeneities of the medium, non linear rheological properties and so on. We reproduce known solutions to test the method and then we furnish new calculations of models which cannot be solved analytically. We consider spherical and non spherical (sill and spheroid) sources. Furthermore, different medi- um properties are included, such as density contrasts, elastic inhomogeneities and rheolog- ical variations. The models are discussed as a preliminary approach to reproduce the gravity changes observed at Campi Flegrei during the 1982-1984 unrest episode. 2. The Campi Flegrei caldera The Campi Flegrei caldera, hereafter CF, near Naples (Italy), is a complex volcanic sys- tem with several, mostly monogenic, explosive vents (fig. 1) in a densely populated area of some 80 km2. The CF is characterized by an outer caldera rim, about 12 km in diameter (open triangles in fig. 1) and an inner caldera of 6-8 km in diameter (full triangles). Slow and re- markable ground movements are typical of this area, as recorded since Roman times. In 1970- 1972 and 1982-1984, a cumulative uplift of about 3.5 m was observed, followed by changes in the shallow hydrothermal system (Chiodini et al., 2003) and an increased rate of shallow seismicity (see De Natale et al., 1991, for an overview). Following January 1985, the ground began to subside at a much slower rate than dur- ing the uplift. The deformation pattern, both during the uplift and the subsidence, exhibits a dominant axial symmetry, concentrated in the inner caldera, with the maximum uplift located at the caldera center (the city of Pozzuoli). Gravity measurements at CF have been carried out since 1981 and the largest gravity variations are also observed at Pozzuoli. The gravity changes are well correlated with the elevations changes (stars) as reported in fig. 2a both dur- ing the 1982-1984 uplift and the 1985-91 defla- Fig. 2a-c. a) gravity variations and elevation change observed during 1981-2001 at Serapeo ro- man market (Pozzuoli), at the center of CF; b) grav- ity due to the ground uplift (free-air correction); c) residual gravity due to the difference between the ob- served and the free-air corrected gravity. tion phase. In particular, during the unrest phase, at Serapeo (a roman market in the center of Pozzuoli) the gravity change normalized to the uplift was −216±7 µGal/m, in good agree- ment with the average of all the stations −213±6 µGal/m. This value is comparable to those measured during inflation episodes at Long Valley and Rabaul calderas. However, during the deflation phase, a ratio between gravity change and uplift of −224±24 µGal/m (similar to that computed during uplift) was observed only at Serapeo (the gravity benchmark closest to the point of maximum uplift), while more Vol51,1,2008_DelNegro 5-02-2009 10:28 Pagina 121 122 E. Trasatti and M. Bonafede scattered ratios were inferred at other stations (e.g., Solfatara and La Pietra), possibly due to local effects (Berrino, 1994). The direct meas- urement of the free-air vertical gradient at CF is γ=−290±5 µGal/m (Berrino et al., 1984), about 20 µGal/m smaller (in absolute value) than the reference value (fig. 2b). The residual gravity ∆gR (the difference between the observed ∆g and the free-air gravity ∆gFA) increased with the uplift by 75±12 µGal/m during 1982-1984 (fig. 2c). The 1982-1984 crises were interpreted as due to an increase in pressure of a deep magma chamber, while geochemical evidence suggested thermodynamic disturbances of shallow aquifers as a response to increased heat flow from below (Bullettin Volcanologique, 1984). The bell shape of the vertical surface displacement suggests modelling the data by means of Mogi isotropic pressure sources. Gottsmann et al. (2006) pre- sented inversions of deformation data adopting different source geometries: point-source, sill and spheroidal cavity. All the sources were locat- ed at about 3 km below sea level, and they sug- gested that the spheroidal source may depict an envelope around a hybrid of both magmatic and hydrothermal sources. Battaglia et al. (2006) performed joint inversions of gravity and defor- mation data, finding that a sill-like source is pre- ferred to interpret the inflation phase, while a shallower spheroid source is more suitable for the deflation phase. The works mentioned above assume the medium to be elastic and homoge- neous. These approximations only provide limit- ed insights into subsurface dynamics at areas such as CF, where intense faulting and inelastic deformations are likely to take place. 3. Method of gravity calculation by FE 3.1. Gravity contributions In general, the total gravity change observed at a benchmark can be separated into the term depending on the elevation change (removed by free-air correction) and the term due to mass re- distribution: (3.1)g gg FA R∆ ∆∆ += where ∆gFA=γuz (γ is the free-air gradient and uz is the vertical displacement of the benchmark). As suggested by many authors, the data is cor- rected using the measured free-air vertical gra- dient (e.g., Berrino, 1994). The residual gravity change is given by (3.2) where G is the gravitational constant, V is the volume over which the density variations ∆ρ do not vanish, rr is the vector between dV and the observation point and θ is the angle between rr and the vertical (fig. 3). According to Sasai (1986) and Bonafede and Mazzanti (1998), ∆ρ in eq. (3.2) accounts for three contributions (3.3) where δρs is the density change related to the introduction of new mass from remote distance, ρ is the material density in the reference config- uration, εij is the strain tensor and u is the dis- placement field. The first term accounts for the new mass intruded; the second term is the rela- tive density change arising from the compress- ibility of the medium; the third one is due to density variations within the medium. Consid- ering the three density changes indicated in eq. (3.3), ∆gR can be written as (fig. 3): (3.4) where ∆gS depends on δρs (vanishing for a massless source) and ∆gM is the contribution due to deformation of the medium surrounding the source. The term ∆gV is expressed by (3.5) and depends on the finite compressibility of rocks; it vanishes if the medium is incompress- ible. The ∆gL term depends on u⋅∇ρ and it ac- counts for the Bouguer correction (expressing the contribution of the lifted mass above the free surface) and the source inflation/deflation (from the displacement of source boundaries). Furthermore, it accounts for the displacement cos g G r dVkk V V 2ρε θ ∆ = − # withg g g g g gR S M M V L∆ ∆ ∆ ∆ ∆ ∆= + = + us kk $4ρ δρ ρε ρ∆ = − − cos g G r dVR V 2ρ θ ∆ ∆= # Vol51,1,2008_DelNegro 5-02-2009 10:29 Pagina 122 123 Gravity changes due to overpressure sources in 3D heterogeneous media: application to Campi Flegrei caldera, Italy of volumes with different density if the medium is not homogeneous. If ρ varies continuously, ∆gL is given by (3.6) If ρ is discontinuous across a surface S within V, the exceeding mass depends on the displacement of the boundary between densities ρin and ρout as (ρout−ρin) u⋅dS. Note that the scalar product be- tween the orientation of the density interface dS and its displacement u accounts for contribu- tions of any orientation of the surface. There- fore, the contribution due to the density discon- tinuities can be expressed as (3.7) where S is the surface over which the density changes. It is important to note that ∆gL may in- clude the contribution of the new mass within the source, if the density gradient is chosen properly. However, we separate the effect of the ( ) cos g G r du Sout inL S 2 $ρ ρ θ ∆ = − # cos g G r dVu V L 2$4ρ θ ∆ = - # displacement of the source boundary (within ∆gL term) from the effect of the new mass (∆gS) because the density contrast between inside/ outside the source is unknown, being one of the main target of gravity modelling. 3.2. Numerical integration of gravity within FE models The subdivision in single contributions of eq. (3.4) is a suitable method to compute gravi- ty changes from FE models of deformation. We consider massless cavities since the source term ∆gS is treated separately in section 5. The grav- ity changes are calculated in two separate steps: 1) FE calculation of displacement u and strain εij fields in the medium, due to the source inflation; 2) computation of gravity terms ∆gL and ∆gV from displacements and strains obtained at 1). Since the gravity terms ∆gL and ∆gV are in- tegrals over the total volume, they are comput- ed as a sum of contributions within each finite element: (3.8) where N is the total number of elements, Vi is the volume of the i-th element of the grid and f is a generic integrand function. The gravity term ∆gV expressed in eq. (3.5) can be easily evaluated within the FE domain as (3.9) where ρi is the density of the i-th element. The term ∆gL due to the density variation is compli- cated by the integration over the faces of ele- ments. Indeed, for each face of the i-th brick el- ement we must consider the density variations within the adjacent element. The term ∆gL from eq. (3.7) is expressed by (3.10) where the surface integral is performed over each j-th face of the i-th element Sij and ρin is the density of the element considered, while ρout is ( ) cos g r du Sout inL ji N S 1 6 1 2 ij $ρ ρ θ ∆ = - == // # cos g G r dVV i i N kk V 1 2 i ρ ε θ ∆ = - = / # ( , , ) ( , , )f x y z dV f x y z dV Vi i N V 1 = = / ## Fig. 3. Schematic representation of the gravity con- tributions: ∆gFA is the free-air change proportional to the uplift; ∆gL is the change caused by the lifted por- tion of the ground (Bouguer anomaly) and of inter- faces between density layers; ∆gV is the gravity field due to the medium compressibility; ∆gS is the contri- bution of the material filling the expanding part of the source. In the right hand side of the figure is in- dicated the geometry of the gravity integration of el- ementary volume dV, its displacement u, the obser- vation point P, the vector R between dV and P, and the angle θ between them. Vol51,1,2008_DelNegro 5-02-2009 10:29 Pagina 123 124 E. Trasatti and M. Bonafede the density of the adjacent element sharing the face Sij. Note that this contribution must be evaluated only once. For the elements forming the free surface, this contribution is the Bouguer correction (since the faces lying on the free surface have confining density equal to 0). If a density contrast is present within two or more elements of the medium, ∆gL accounts for this contribution. The great advantage of this integration is that non planar interfaces between layers can be considered. Each element may have arbitrarily shape and may be characterized by different density, elastic parameters or rheol- ogy independently from the rest of the model. The integrals over single elements are eval- uated numerically using the Gauss quadrature formula. We give a brief overview of this tech- nique: 1) the Gauss integration is computed by mapping each arbitrarily distorted brick ele- ment, the (x, y, z) space, to a trilinear hexahe- dral element of side 2l, the (xl, yl, zl) space. The transformation applies to local cells only (fig. 4a). The volume integration for an element be- comes (3.11) where J is the Jacobian determinant. 2) The volume integral is evaluated as the sum of the integrand function calculated at Gauss points. We adopt the two-points integra- tion (for each dimension), which requires to ap- proximate the integrand function at 8 Gauss points: .( , , ) , ,x y z l l l3 3 3G G G ! ! !=l l l _ i ( , , ),y x y z f( , , ) ( , , ), ( , , ) f x y z dV x x y z z x y z V l l l l l l = −−− ( , , )J x y z dV l l l l l l l l l l l l l 6 @ # ### Fig. 4a,b. a) the Gauss integration is computed by mapping each arbitrarily distorted brick element (x-space) to a trilinear hexahedral element (xl-space). The two-points Gauss quadrature is performed by calculating the in- tegrand function at the 8 Gauss points indicated by «+» symbol. b) The surface integral within element faces is performed by mapping x into xl. The vectors ds1 and ds2 are defined to be tangent to the surface, defining the unit area shaded (see text for details). Vol51,1,2008_DelNegro 5-02-2009 10:29 Pagina 124 125 Gravity changes due to overpressure sources in 3D heterogeneous media: application to Campi Flegrei caldera, Italy The volume integral becomes the sum of the fol- lowing values: (3.12) Since the Gauss integration is the technique adopted by FE to calculate many derived fields such as strain and stress, these values are direct- ly available from our code (MARC, 1994, ver- sion 2005) at Gauss points without any further approximations. The term ∆gV in eq. (3.9) can be computed from eq. (3.12). However, compu- tation of ∆gL in eq. (3.10) involves surface inte- grals over each element face. Following the de- velopment outlined in Greenberg (1978), each face is locally transformed in space (xl, yl, zl), defining the vectors ds1 and ds2 tangent to the plane considered. In the example shown in fig. 4b, the face 4-3-7-8 is mapped into the xl-con- stant plane. The vector ds1 is defined to be along zl-constant curve and ds2 is along yl-con- stant curve. The elemental area vector dS, de- noted by the shaded parallelogram, can be com- puted as the vector product between ds2 and ds1 (positive pointing out of the cell). The surface integral becomes (13) and it is computed by applying the two-points Gauss rule in two-dimensions. Indeed, 4 Gauss points are determined over the mapped plane, e.g., , and the integrand function is approximated at these points follow- ing the procedure outlined above. 4. Cases of study We develop an FE model characterized by about 75000 8-nodes brick elements, extending 160×160×80 km3 (fig. 5). The domain is large enough to avoid bias from boundaries over which displacement and stress fields are im- posed to vanish. Furthermore, the bottom boundary is kept fixed. The grid is character- ized by a flat interface at 3 km depth which may be used to represent a shallow layer, and a , ,l l l3 3! !_ i ( )u S u s sf d f d d1 2 S S $ $ #=# # f ( , , ) ( , , )x y z f x y z l l l l l l G G G 1 8 = --- dVl l l l l l l/### cylindrical boundary of radius 3 km and depth 3 km lying at the center of the model, to ap- proximate the CF caldera. The element size is about 350 m side in the central part of the grid, while it increases with distance up to 10 km close to the domain boundaries. We show examples of gravity computations for pressurized sources of various shapes in ho- mogeneous/inhomogeneous media. In all the cases shown in this section, the sources are con- sidered empty so that only ∆gM is computed. A brief outline of the main characteristics of the models are reported in table I. Unless changes are due to specific model configurations, the medium is considered homogeneous, elastic and isotropic with rigidity µ0=1 GPa and Pois- son ratio ν= =0.25. All the sources are placed at 5 km depth, at the center of the model, co-axial with the caldera and undergo to overpressure ∆P=50 MPa. When the density is homoge- neous, its value is ρ0=2500 kg/m3. The models named MOGI-n are characterized by a spheri- cal source of radius a=1 km; in particular the case with homogeneous and elastic medium is MOGI-1. Density contrasts are considered in models MOGI-2 and MOGI-3, characterized by Fig. 5. Perspective view of the FE model. The caldera is a cylinder of radius 3 km and height 3 km, surrounded by a layer of similar thickness on top of the halfspace. The Mogi source is shown at a depth of 5 km, with a radius of 1 km. Vol51,1,2008_DelNegro 5-02-2009 10:29 Pagina 125 126 E. Trasatti and M. Bonafede a density equal to ρ1=1800 kg/m3 within the shallow layer and the caldera (respectively); the density of the remaining medium is ρ0. In mod- el MOGI-4 the caldera has heterogeneous rigid- ity with respect to the remaining medium µ1= =0.1µ0, while the density ρ0 is constant any- where. Model MOGI-5 contains both the last characteristics: the caldera has heterogeneous rigidity µ1 and density ρ1, while the rest of the medium has rigidity µ0 and density ρ0. We also consider two inelastic models, assuming a plas- tic rheology of the medium to describe the highly fractured rocks of the CF area. The medium behaves plastically when the maxi- Table I. List of models of which we calculate the gravity variations. The overpressure within the sources is ∆P=50 MPa and all the sources are placed at depth d = 5 km. The layer has horizontal interface and a height of 3 km. The caldera has a cylindrical shape of 3 km radius and 3 km height; it lays at the center of the model, co- axial with the sources. The parameters are explained in the text. Model Source Medium characteristics Parameters MOGI-1 Sphere Homogeneous and elastic µ0 = 1 GPa; ρ0 = 2500 kg/m3; a = 1 km MOGI-2 ” Heterogeneous density layer ρ1 = 1800 kg/m3; ρ0 = 2500 kg/m3 MOGI-3 ” Heterogeneous density caldera ” MOGI-4 ” Heterogeneous rigidity caldera µ1 = 0.1 µ0 MOGI-5 ” Heterog. rigidity and density cal. µ1 =0.1 µ0; ρ1 =1800kg/m3; ρ0=2500kg/m3 MOGI-6 ” All plastic σy = 30 MPa MOGI-7 ” Plastic caldera σy = 0.5 MPa SPHEROID Spheroid Homogeneous and elastic Vsphere = Vspheroid; a=1850 m; b=c=0.4a SILL Sill ” l = 2 km Fig. 6. Contributions and medium gravity change due to an inflating spherical pressure source in a homoge- neous and elastic half-space. Comparison between analytical solutions (solid line) and numerical (squares). The contribution due to the inflation of a massless sphere ∆gM = ∆gL + ∆gV is null. Vol51,1,2008_DelNegro 5-02-2009 10:29 Pagina 126 127 Gravity changes due to overpressure sources in 3D heterogeneous media: application to Campi Flegrei caldera, Italy mum shear stress is greater than a specified yield stress. The plastic rheology is only treat- able by means of numerical methods, and easi- ly implementable in FE models. In model MO- GI-6 the whole medium behaves plastically with yield stress σy=30 MPa, while MOGI-7 is characterized by plastic caldera (σy=0.5 MPa) and remaining elastic medium, to evidence the highly deformable inner caldera. Two further models are considered, with different shape of the source: spheroid and sill. The spheroid has the same volume of the sphere and aspect ratio equal to 0.4 with resulting semi-major axis a = 1850 m and semi-minor axes b = c = 740 m; it is vertically elongated. The sill is a horizontal square crack of side l = 2 km; both sources act in a homogeneous and elastic medium, to consider the effects of the different shape only. Figure 6 reports the gravity patterns from the source axis to 10 km distance along the x-axis; the solid line is the analytical solution while the squares are the results from numerical integra- tion. The figure shows that the gravity change ∆gM due to a Mogi source in a homogeneous and elastic half-space is null and the only contribu- tion observable at surface would be the free-air effect. This agrees with analytical results by Walsh and Rice (1979), Sasai (1986), Bonafede and Mazzanti (1998). Our numerical «zero» is 0.5 µGal, well below instrumental error, con- firming the robustness of the Gauss integration within the FE mesh and the good approximation of the half-space by the numerical domain. Different results are obtained for models SPHEROID and SILL, dilating in a homoge- neous medium (fig. 7a,b). We report the comput- ed vertical displacement (the thin solid line) which is representative of ∆gFA. From now on the gravity changes will be shown normalized to the uplift above the source center (coinciding with the maximum up lift for sphere and sill models). In the case of the spheroid, the gravity variation ∆gM is negative (fig. 7a, thick solid line). Therefore, if the gravity signal generated by a spheroid is modeled by a spherical source, it may lead to an overestimation of the density within the source. This result is also confirmed by Battaglia and Segall (2004) that performed in- versions of geodetic and gravity data generated by a spheroid, by means of a spherical source. They found that the isotropic point-source ap- proximation leads to an overestimation of depth, mass and density of the intrusion. Fig. 7b shows the opposite effect due to the sill: the total contri- bution is positive, without new mass input. Note that in this case the contribution ∆gV (dotted line) is almost null. This behavior arises from the strain field generated by the sill in the half-space, causing contraction above it and lateral dilata- tion. We compare these results with the analyti- cal ones provided by Sasai (1986), finding good agreement (G. Currenti personal communica- Fig. 7a,b. Gravity terms normalized to the maximum uplift uzmax, and vertical displacement in a homogeneous medium due to (a) vertically elongated spheroidal source, (b) horizontal sill (in this case ∆gV 0 and ∆gM over- laps ∆gL). - Vol51,1,2008_DelNegro 5-02-2009 10:29 Pagina 127 128 E. Trasatti and M. Bonafede tion). In a recent paper Battaglia et al. (2006) im- plemented the analytical solutions of the sphe- roid (Yang et al., 1988) within the gravity calcu- lation method developed by Walsh and Rice (1979). They found that the gravity change due to the medium deformation is negligible for a spheroidal cavity while it is noticeable for a sill. Good agreement is found from the comparison between our and their results for the sill. Instead, our results for the spheroidal cavity is not in agreement with their findings, as shown in fig. 7a. However, since comparisons between our models and the analytical solutions (when avail- able, like models MOGI-1, MOGI-2 and SILL) show very good agreement, and since the method of gravity computation shown here is in- dependent on the source shape and properties of the medium, we are confident with the robust- ness of the present results. Figure 8a,b reports the results for two models characterized respectively by density contrast of 30% at 3 km depth (MOGI-2) and within the Fig. 8a,b. Gravity terms normalized to the maximum uplift uzmax, and vertical displacement due to an inflating spherical pressure source in a medium characterized by density contrasts of 30% in: a) shallow layer of height 3 km; b) cylindrical caldera. Fig. 9a,b. Gravity terms normalized to the maximum uplift uzmax, and vertical displacement due to an inflating spherical pressure source in a medium characterized by: a) heterogeneous caldera µ1/µ0 = 0.1; b) caldera with heterogeneous rigidity µ1 and density ρ1. Vol51,1,2008_DelNegro 5-02-2009 10:29 Pagina 128 129 Gravity changes due to overpressure sources in 3D heterogeneous media: application to Campi Flegrei caldera, Italy caldera rim (MOGI-3). In the first case (fig. 8a), solutions are in agreement with those by Bonafede and Mazzanti (1998) (not shown here). A larger effect is visible when the low density is restricted to the caldera (fig. 8b). In this case, ∆gL accounts for the density discontinuity at the cir- cular boundary of the cylinder also. The density jump reflects into the step of the ∆gL pattern, while the volumetric and free-air (∝ uz) terms are continuous. Indeed, if we compare the two mod- els, we find a similar pattern of ∆gV because the density contrast has little effect in the strain field, while ∆gL is more sensitive to this discontinuity. These models show that this kind of inhomo- geneities, even if very localized (like within the caldera), must be taken into account to perform gravity modelling. Results similar to MOGI-2 were obtained by Battaglia and Segall (2004) for a massless cavity in a layered medium. They found a small (but different from 0) contribution to the residual gravity due to a density contrast of 13% at Long Valley caldera. The contribution of heterogeneous rigidities within the medium is emphasized in model MO- GI-4. Indeed, as long as the rigidity ratios are small, the contribution to ∆gM/uz is of the order of some µGal/m. In fig. 9a the results is shown for an extreme ratio of µ1/µ0=0.1 between caldera and remaining medium. The low rigidity of the caldera amplifies the displacement field and hence the terms linked to this parameter (∆gFA and ∆gL), while the volumetric term ∆gV is al- most unchanged from the homogeneous case in fig. 6. A very low effective rigidity could be suit- able to describe longterm deformation for a stan- dard linear solid rheology; this rheology may ap- ply to volcanic areas such as CF, due to the high temperature, erupted products, incoherent materi- als and hydrothermal activity. A further model, MOGI-5, accounts for both the characteristics of MOGI-4 and MOGI-3, having the caldera with low rigidity and low density (fig. 9b). Since the displacement generated from this model is equal to MOGI-4, ∆gV is very similar to that of panel (a). The density contrast between the caldera and the remaining medium causes a lower ∆gL and a small step in its pattern. We also consider two cases of inelastic medium, described by plastic rheology. In the case of uniform yield stress within the medium, the deformation is approximately radial around the spherical source and the total effect is null (fig. 10a), similarly to the elastic homogeneous model. The sill is also modelled (not shown here) in a medium homogeneously plastic, find- ing that ∆gM/uz is unchanged with respect to the elastic case, as observed for the spherical source. Indeed, the plastic medium contributes by amplifying both uz and the gravity terms, leaving the ratio between them unchanged. In the case of model MOGI-7, a rheologic discontinuity is present, since the caldera is plas- tic while the rest of the medium is elastic. It is ev- ident from fig. 10b that when the yield stress is not uniform, the gravity change ∆gM is different from 0. Indeed, the total contribution is positive within the caldera while it becomes negative out of the caldera rim, vanishing at 7-8 km from the center of the model. This model is useful to de- scribe localized rheology discontinuities, likely to be present in volcanic areas. It is interesting to note that the gravity terms are very similar to those shown in fig. 9a, indicating that both the low rigidity and the plastic rheology act to en- hance the deformation within the caldera, allow- ing for larger ∆gL terms while ∆gV is similar to the elastic value. Indeed, ∆gL is sensitive to the deviatoric strain (mostly plastic) while ∆gV is sensitive only to the isotropic strain (elastic). This kind of behavior was also figured out by Trasatti et al. (2005), showing that the elasto-plastic rhe- ology induces larger deformations within the caldera. 5. Contribution of the inflating source The previous section considered the reser- voirs massless, i.e. cavities within which the overpressure ∆P is assigned. Whether the ex- pansion of cavities is due to new mass input, in- ternal processes or pore pressure and tempera- ture migration, the gravity change due to the source at a benchmark can be computed as (5.1) cos cos g G r dV r dV S S V S V 2 2 S S ρ θ ρ θ ∆ = - -l l l l l = G # # Vol51,1,2008_DelNegro 5-02-2009 10:29 Pagina 129 130 E. Trasatti and M. Bonafede where VlS, VS are the source volumes after and before the deformation, respectively, and all primed variables refer to the deformed configu- ration. We assume that the main process deter- mining the ∆gS is the new mass entering into the reservoir. This process has several conse- quences: volume variation, density growth due the increased internal pressure and the displace- ment of the center of mass of the reservoir due to the asymmetric deformation of the medium in proximity of the free surface. For simplicity, we assume ρS, ρlS to be mean internal densities, and we approximate the mass at the source cen- ter (due to the source depth, assumed as 5 km, much larger than the extension). The mass change is ∆M=ρ∆V +V∆ρ and the density after deformation is expressed by Fig. 10a,b. Gravity terms normalized to the maximum uplift uzmax, and vertical displacement in an inelastic medium: a) all plastic, b) plastic caldera and remaining elastic. Fig. 11. Density calculations versus volume variations of the proposed models. Vol51,1,2008_DelNegro 5-02-2009 10:29 Pagina 130 131 Gravity changes due to overpressure sources in 3D heterogeneous media: application to Campi Flegrei caldera, Italy ρlS=(M + ∆M)/(V + ∆M), at the first order. It re- sults that (5.2) The ratio ∆ρ/ρS=∆P/K∼10−3, where ∆P is the overpressure and K is the isothermal incom- pressibility, so that ρlS ∼ ρS and the main effect on the gravity change is the volume increase of the magma reservoir (Franchini, 2005). In con- clusion, the resulting ∆gS is expressed by (5.3) The volume variation of a cavity with arbitrary shape is obtained by integrating the normal dis- placement over the boundary (5.4) Within the FE grid the integral is performed adopting the procedure outlined for the gravity term ∆gL (see eqs. (3.7) and (3.10)), using eq. (3.13). The observed gravity change at CF is 75±12 µGal/m, which must be equal to ∆gR= =∆gS + ∆gM. From eq. (5.3) and the values of ∆gM computed for the sources in table I we ob- tain the source densities ρS. Volume variations and density estimates are reported in fig. 11. The density value calculated V du S S $∆ = # cos g G V r s S 2ρ θ ∆ ∆= 1 V M V V M V V V V V 1S S S S + +ρ ρ ρ ρ ρ ρ ∆ ∆ ∆ ∆ ∆ ∆ − − + + + = + +l b bl l for the sphere in the homogeneous medium (MOGI-1) is very similar to that usually found in literature for magmatic sources. Indeed, Berrino et al. (1984) and Berrino (1994) esti- mated 2500 kg/m3, suggesting the hypothesis of silicate melts entering a reservoir rather than hydrothermal fluids within a confined aquifer. All the models (except MOGI-6) show similar volume variation, while the densities inferred range from ∼400 to ∼3300 kg/m3. The contribu- tion of the medium, when homogeneously plas- tic (MOGI-6), is null and the inferred density is comparable to the elastic models, ∼2000 kg/m3. However, the volume variation is very large, due to the enhanced deformation allowed in the medium. The density inferred in models char- acterized by density layering (MOGI-2 and MOGI-3) is greater than that inferred in a ho- mogeneous model. Indeed, since the gravity contribution ∆gM is negative, the observed pos- itive gravity change must be due to a larger den- sity for the new mass. Also MOGI-5 (low rigid- ity and low density caldera) requires high den- sity in spite of the positive ∆gM value, because of the larger uplift. It is interesting to note that models MOGI-4 (low rigidity caldera), MOGI- 5 and MOGI-7 (plastic caldera) produce very similar ∆gM value (fig. 9 and 10), but intrusion density differ by ±25%. Due to the compensation of the negative value computed for ∆gM (fig. 7a), the density inferred for the spheroid is quite high, being about 3300 kg/m3. This value is much higher Table II. Volume variation and density calculated for each model (characteristics reported in table I) due to in- put of new mass within the cavity. Model ∆V (106 m3) ρS (kg/m3) ∆M (1011 kg) MOGI-1 155 2690 4.2 MOGI-2 155 2907 3.8 MOGI-3 155 3194 3.4 MOGI-4 159 2523 4.0 MOGI-5 159 3328 5.3 MOGI-6 354 2033 7.2 MOGI-7 159 2031 3.2 SPHEROID 174 3271 5.7 SILL 125 410 0.5 Vol51,1,2008_DelNegro 5-02-2009 10:29 Pagina 131 132 E. Trasatti and M. Bonafede than pertinent to trachytic basalts at CF (~2300- 2500 kg/m3). On the other hand, the density in- ferred for the sill model is very low, 410 kg/m3. Since supercritical fluid densities are of the or- der of hundreds of kg/m3, this value would indi- cate a hydrothermal source instead of a mag- matic source, as suggested by Battaglia et al. (2006) for a penny shaped crack. Computations for mass input ∆M=ρ∆V are reported in table II. Based on the evaluation of gravity residuals and the assumption of a pres- surized point-source, Berrino et al. (1984) and Berrino (1994) suggested that a subsurface mass increase of 2 · 1011 kg took place between 1982-1984. Results from models with spherical sources are typically 2-3 times greater than this value. The all plastic model, MOGI-6, has a very high volume variation due to the enhanced deformation of the plastic medium. The com- puted density is ∼2000 kg/m3 while the mass in- put is 3-4 times the estimated value. On the oth- er side, the mass calculated for the sill model is a factor 4 less than the reference value. 6. Conclusions We develop a numerical technique to calcu- late gravity changes due to deep inflating sources using FE modelling. Even if the models presented are too simple for realistic application to the CF caldera, the method is very general and accounts for many characteristics of the studied area: various source shapes, density and rigidity discontinuities, rheological hetero- geneities. Very few analytical solutions are available in literature (e.g., sphere, sill, spheroid in elastic media and point-sources in viscoelas- tic media), while the FE models may be charac- terized arbitrarily. The gravity calculations dis- tinguish between the contribution of the medi- um and of the mass intrusion. Very different contributions come from the deformation of the medium, according to source shape. We sepa- rate the effect of the displacement of layers in- terfaces (∆gL) and of the medium compressibil- ity (∆gV). We consider the gravity changes oc- curred at CF caldera during the 1982-1984 un- rest episode, where a free-air corrected gravity residual of 75±12 µGal/m was observed. The FE models developed includes heterogeneities in density, rigidity and rheological properties in order to point out which may be important in af- fecting the observations. It must be noted that the large deformation is well outside the elastic limit of any rock material. We show that non- spherical sources such as sill and spheroid yield positive or negative gravity changes (respective- ly), without input of new mass. The density val- ues computed for these models are very differ- ent, changing by a factor 8 between extreme end models. The results clearly point out the impor- tance of the source shape and of the characteris- tics of the medium in gravity calculation and source density estimation. The proposed models, even if not yet orient- ed yet to data inversions, give hints in under- standing the effects of considering the realistic structure of the studied area in gravity computa- tions. We fit the uplift and the gravity above the source but do not yet attempt to reproduce the spatial pattern away from the axis. This will be done after developing a realistic model taking into account the characteristics of the area. Den- sity and elastic parameters can be extracted, in principle, from seismic tomography studies and from empirical relationships between seismic velocity and density (Brocher , 2005). However, long term deformation in volcanic areas may be largely affected by anelastic relaxation and drainage conditions. The large deformations which took place at CF are better described in terms of plastic rather than elastic constitutive relationships. Laboratory studies on samples ex- tracted from deep drilling seem the only plausi- ble way to obtain information on the rheological properties of the medium. The source shape is also very important, but only direct modelling has been attempted up to now. An inversion pro- cedure to retrieve the shape of a point-like de- formation source in laterally heterogeneous and inelastic media have been proposed (Trasatti et al., 2008), and this technique should be further developed to include modeling of gravity changes. From direct modeling we know that vertical and horizontal deformation patterns are necessary to distinguish among different source shapes (e.g., Dieterich and Decker, 1975): in- cluding gravity data gives further constraints to restrict the range of acceptable models. Vol51,1,2008_DelNegro 5-02-2009 10:29 Pagina 132 133 Gravity changes due to overpressure sources in 3D heterogeneous media: application to Campi Flegrei caldera, Italy Acknowledgements This work was supported by funding from Italian Civil Protection INGV-DPC V3-2 Campi Flegrei Project.We thank C. Giunchi for the useful discussions and hints and G. Berrino for providing deformation and gravity data. Useful suggestions from Frank Roth and Gilda Currenti are greatly acknowledged. REFERENCES BATTAGLIA, M., and P. SEGALL (2004): The interpretation of gravity changes and crustal deformation in active vol- canic areas, Pure Appl. Geophys., doi: 10.1007/ s00024-004-2514-5. BATTAGLIA, M., C. TROISE, F. OBRIZZO, F. PINGUE and G. De NATALE (2006): Evidence for fluid migration as the source of deformation at Campi Flegrei caldera (Italy), Geophys. Res. Lett., doi: 10.1029/2005GL024904. BERRINO, G. (1994): Gravity changes induced by height- mass variations at Campi Flegrei caldera, J. Volcanol. Geotherm. Res., 61, 293-309. BERRINO, G., G. CORRADO, G. LUONGO and B. TORO (1984): Ground deformation and gravity change accompanying the 1984 Pozzuoli uplift, Bull. Volcanol., 47, 187-200. BONAFEDE, M. and M. MAZZANTI (1998): Modelling gravi- ty variations consistent with ground deformation in the Campi Flegrei caldera (Italy), J. Volcanol. Geotherm. Res., 81, 137-157. BROCHER, T.M. (2005): Empirical relations between elastic wavespeeds and density in the earth’s crust, Bull. Seism. Soc. Am., 95, 2081-2092. BULLETTIN VOLCANOLOGIQUE (1984): The 1982-1984 bradyseismic crisis at Campi Flegrei, Italy, Bull. Vol- canol., 47. CHIODINI, G., M. TODESCO, S. CALIRO, C. DEL GAUDIO, G. MACEDONIO and M. RUSSO (2003): Magma degassing as a trigger of bradyseismic event: the case of Phlegre- an Fields (Italy), Geophys. Res. Lett., 30, doi: 10.1029/ 2002GL016790. CURRENTI, G., C. DEL NEGRO and G. GANCI (2007): Model- ling of ground deformation and gravity fields using fi- nite element method: an application to Etna volcano, Geophys. J. Int., 169, 775-786, doi: 10.1111/j.1365- 246X.2007.03380.x. DE NATALE, G., F. PINGUE, P. ALLARD and A. ZOLLO (1991): Geophysical and geochemical modelling of the 1982- 1984 unrest phenomena at Campi Flegrei caldera (South- ern Italy), J. Volcanol. Geotherm. Res., 48, 199-222. DIETERICH, J.H. and R.W. DECKER (1975): Finite element modelling of surface deformation associated with vol- canism, J. Geophys. Res., 80, 4094-4102. FERNANDEZ, J. and J.B. RUNDLE (1994): Gravity changes and deformation due to a magmatic intrusion in a two lay- ered crustal model, J. Geophys. Res., 99, 2737-2746. FERNANDEZ, J., K.F. TIAMPO and J.B. RUNDLE (2001): Vis- coelastic displacement and gravity changes due to point magmatic intrusions in a gravitational layered solid Earth, Geophys. J. Int., 146, 155-170. FRANCHINI, L. (2005): Deformazioni e variazioni della gravità in aree vulcaniche: influenza della struttura re- ologica, in italian, Università degli Studi di Bologna, Corso di Laurea Specialistica in Fisica. GOTTSMANN, J., H. RYMER and G. BERRINO (2006): Unrest at the Campi Flegrei caldera (Italy): a critical evalua- tion of source parameters from geodetic data inversion, J. Volcanol. Geotherm. Res., 150, 132-145. GREENBERG, M.D. (1978): Foundations of Applied Mathe- matics (Prentice-Hall, Inc., Englewood Cliffs, New Jer- sey), pp. 456. MARC (1994): Analysis Research Corporation (Palo Alto, CA). SASAI, Y. (1986): Multiple tension-crack model for dilatan- cy: surface displacement, gravity and magnetic change, Bull. Earth. Res. Inst., 61, 429-473. TRASATTI, E., C. GIUNCHI and M. BONAFEDE (2005): Struc- tural and rheological constraints on source depth and overpressure estimates at the Campi Flegrei caldera, Italy, J. Volcanol. Geotherm. Res., 144, 105-118. TRASATTI, E., C. GIUNCHI and N.P. AGOSTINETTI (2008): Nu- merical inversion of deformation caused by pressure sources: application to Mt. Etna (Italy), Geophys. J. Int., 172, 873-884, doi: 10.111/j.1365-246x.2007.03677.x. WALSH, J.B. and J.R. RICE (1979): Local changes in gravi- ty resulting from deformation, J. Geophys. Res., 84, 165-170. WANG, R.F., F.L. MARTÌN and F. ROTH (2006), PSGRN/ PSCMP a new code for calculating co- and post-seis- mic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation theory, Com- put. Geosci., 32, 527-541. YANG, X., P.M. DAVIS and J.H. DIETERICH (1988): Deforma- tion from inflation of a dipping finite prolate spheroid in an elastic halfspace as a model for volcanic stress- ing, J. Geophys. Res., 93, 4249-4257. Vol51,1,2008_DelNegro 5-02-2009 10:34 Pagina 133 Vol51,1,2008_DelNegro 5-02-2009 10:34 Pagina 134