Gravity changes due to overpressure sources in 3 D heterogeneous media : application to Campi Flegrei caldera , Italy

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 media. 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 contributions also come from density and rigidity heterogeneities within the medium. Furthermore, if the caldera is elasto-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/m (super critical water) to ∼3300 kg/m (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


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 mapping subsurface density distributions.Both of them contribute to quantify the change in subsurface mass.Several models have been developed 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 redistribution 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) suggested the separation of density variations affecting 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 (1994) developed a method to compute gravity variations due to point magma intrusion in a horizontally layered elastic-gravitational media.The solutions were extended and generalized to point-sources in multi-layered viscoelastic-gravitational media (Fernandez et al., 2001).A comprehensive tool for evaluating deformation and gravity changes due to dislocation in horizontally layered viscoelastic media was also provided recently by Wang et al. (2006).The approaches described above are only partially able to simulate the complex geology of volcanic regions, since they model pointsources in horizontally layered media.Furthermore, 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 maximum values within the caldera.This behavior allows us to estimate a source depth of 5 km (plastic medium) instead of 3 km (elastic medium), adopting the same overpressure and in accordance with information from groundwater studies and petrological data.Recently, Currenti 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 neglected.
In the present paper we attempt to model both the deformation and the gravity.We develop a technique to calculate gravity variations due to deep pressurized sources in heterogeneous media.The method consists in two separated steps: i) the displacement and strain fields due to the source are computed by FE model- 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: arbitrary 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 medium properties are included, such as density contrasts, elastic inhomogeneities and rheological variations.The models are discussed as a preliminary approach to reproduce the gravity changes observed at Campi Flegrei during the 1982-1984 unrest episode.

The Campi Flegrei caldera
The Campi Flegrei caldera, hereafter CF, near Naples (Italy), is a complex volcanic system with several, mostly monogenic, explosive vents (fig. 1) in a densely populated area of some 80 km 2 .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 remarkable ground movements are typical of this area, as recorded since Roman times.In 1970In -1972In and 1982In -1984, a , 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 during 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  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 agreement 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 scattered ratios were inferred at other stations (e.g., Solfatara and La Pietra), possibly due to local effects (Berrino, 1994).The direct measurement 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) presented inversions of deformation data adopting different source geometries: point-source, sill and spheroidal cavity.All the sources were located at about 3 km below sea level, and they suggested 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 deformation data, finding that a sill-like source is preferred 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 homogeneous.These approximations only provide limited insights into subsurface dynamics at areas such as CF, where intense faulting and inelastic deformations are likely to take place.

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 redistribution: 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 corrected using the measured free-air vertical gradient (e.g., Berrino, 1994).The residual gravity change is given by where G is the gravitational constant, V is the volume over which the density variations ∆ρ do not vanish, r r is the vector between dV and the observation point and θ is the angle between r r and the vertical (fig.3).According to Sasai (1986) and Bonafede and Mazzanti (1998), ∆ρ in eq.(3.2) accounts for three contributions where δρs is the density change related to the introduction of new mass from remote distance, ρ is the material density in the reference configuration, εij is the strain tensor and u is the displacement field.The first term accounts for the new mass intruded; the second term is the relative density change arising from the compressibility of the medium; the third one is due to density variations within the medium.Considering 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 incompressible.The ∆gL term depends on u⋅∇ρ and it accounts 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 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 between the orientation of the density interface dS and its displacement u accounts for contributions of any orientation of the surface.Therefore, the contribution due to the density discontinuities can be expressed as where S is the surface over which the density changes.It is important to note that ∆gL may include the contribution of the new mass within the source, if the density gradient is chosen properly.However, we separate the effect of the 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.

Numerical integration of gravity within FE models
The subdivision in single contributions of eq. ( 3.4) is a suitable method to compute gravity changes from FE models of deformation.We consider massless cavities since the source term ∆g S is treated separately in section 5.The gravity 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 ∆g L and ∆gV from displacements and strains obtained at 1).
Since the gravity terms ∆gL and ∆gV are integrals over the total volume, they are computed as a sum of contributions within each finite element: where N is the total number of elements, V i 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 complicated by the integration over the faces of elements.Indeed, for each face of the i-th brick element we must consider the density variations within the adjacent element.
The term ∆gL from eq. (3.7) is expressed by 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 3. Schematic representation of the gravity contributions: ∆gFA is the free-air change proportional to the uplift; ∆gL is the change caused by the lifted portion of the ground (Bouguer anomaly) and of interfaces between density layers; ∆gV is the gravity field due to the medium compressibility; ∆gS is the contribution of the material filling the expanding part of the source.In the right hand side of the figure is indicated the geometry of the gravity integration of elementary volume dV, its displacement u, the observation point P, the vector R between dV and P, and the angle θ between them.
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 rheology independently from the rest of the model.
The integrals over single elements are evaluated numerically using the Gauss quadrature formula.We give a brief overview of this technique: 1) the Gauss integration is computed by mapping each arbitrarily distorted brick element, the (x, y, z) space, to a trilinear hexahedral element of side 2l, the (xl, yl, zl) space.The transformation applies to local cells only (fig.4a).The volume integration for an element becomes 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 integration (for each dimension), which requires to approximate the integrand function at 8 Gauss points: . The volume integral becomes the sum of the following 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 directly available from our code (MARC, 1994(MARC, , version 2005) ) at Gauss points without any further approximations.The term ∆gV in eq.(3.9) can be computed from eq. (3.12).However, computation of ∆g L in eq.(3.10) involves surface integrals over each element face.Following the development outlined in Greenberg (1978), each face is locally transformed in space (xl, yl, zl), defining the vectors ds 1 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-constant plane.The vector ds1 is defined to be along zl-constant curve and ds2 is along yl-constant curve.The elemental area vector dS, denoted by the shaded parallelogram, can be computed as the vector product between ds 2 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 following the procedure outlined above.

Cases of study
We develop an FE model characterized by about 75000 8-nodes brick elements, extending 160×160×80 km 3 (fig.5).The domain is large enough to avoid bias from boundaries over which displacement and stress fields are imposed to vanish.Furthermore, the bottom boundary is kept fixed.The grid is characterized by a flat interface at 3 km depth which may be used to represent a shallow layer, and a cylindrical boundary of radius 3 km and depth 3 km lying at the center of the model, to approximate 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 homogeneous/inhomogeneous media.In all the cases shown in this section, the sources are considered 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 Poisson 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 homogeneous, its value is ρ 0=2500 kg/m 3 .The models named MOGI-n are characterized by a spherical 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 a density equal to ρ1=1800 kg/m 3 within the shallow layer and the caldera (respectively); the density of the remaining medium is ρ0.In model MOGI-4 the caldera has heterogeneous rigidity with respect to the remaining medium µ1= =0.1µ0, while the density ρ0 is constant anywhere.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 plastic 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, coaxial with the sources.The parameters are explained in the text.mum shear stress is greater than a specified yield stress.The plastic rheology is only treatable by means of numerical methods, and easily 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.

Model
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 integration.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 contribution 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, confirming 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 homogeneous medium (fig.7a,b).We report the computed 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 inversions of geodetic and gravity data generated by a spheroid, by means of a spherical source.They found that the isotropic point-source approximation leads to an overestimation of depth, mass and density of the intrusion.Fig. 7b shows the opposite effect due to the sill: the total contribution is positive, without new mass input.Note that in this case the contribution ∆g V (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 dilatation.We compare these results with the analytical ones provided by Sasai (1986), finding good agreement (G.Currenti personal communica-  tion).In a recent paper Battaglia et al. (2006) implemented the analytical solutions of the spheroid (Yang et al., 1988) within the gravity calculation 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 available, like models MOGI-1, MOGI-2 and SILL) show very good agreement, and since the method of gravity computation shown here is independent on the source shape and properties of the medium, we are confident with the robustness 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 uz max , 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 uz max , 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.
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, ∆g L accounts for the density discontinuity at the circular 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 models, we find a similar pattern of ∆gV because the density contrast has little effect in the strain field, while ∆g L is more sensitive to this discontinuity.These models show that this kind of inhomogeneities, 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 ∆g M/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 almost unchanged from the homogeneous case in fig.6.A very low effective rigidity could be suitable to describe longterm deformation for a standard linear solid rheology; this rheology may apply to volcanic areas such as CF, due to the high temperature, erupted products, incoherent materials 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, ∆g V 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, finding that ∆g M/uz is unchanged with respect to the elastic case, as observed for the spherical source.Indeed, the plastic medium contributes by amplifying both u z 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 plastic while the rest of the medium is elastic.It is evident 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 describe 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 enhance the deformation within the caldera, allowing 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 rheology induces larger deformations within the caldera.

Contribution of the inflating source
The previous section considered the reservoirs massless, i.e. cavities within which the overpressure ∆P is assigned.Whether the expansion of cavities is due to new mass input, internal processes or pore pressure and temperature migration, the gravity change due to the source at a benchmark can be computed as where Vl S, VS are the source volumes after and before the deformation, respectively, and all primed variables refer to the deformed configuration.We assume that the main process determining the ∆gS is the new mass entering into the reservoir.This process has several consequences: 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, ρl S to be mean internal densities, and we approximate the mass at the source center (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  The ratio ∆ρ/ρ S=∆P/K∼10 −3 , where ∆P is the overpressure and K is the isothermal incompressibility, so that ρl S ∼ ρS and the main effect on the gravity change is the volume increase of the magma reservoir (Franchini, 2005).In conclusion, the resulting ∆g S is expressed by The volume variation of a cavity with arbitrary shape is obtained by integrating the normal displacement 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 obtain the source densities ρS.
Volume variations and density estimates are reported in fig.11.The density value calculated 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) estimated 2500 kg/m 3 , 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/m 3 .The contribution of the medium, when homogeneously plastic (MOGI-6), is null and the inferred density is comparable to the elastic models, ∼2000 kg/m 3 .However, the volume variation is very large, due to the enhanced deformation allowed in the medium.The density inferred in models characterized by density layering (MOGI-2 and MOGI-3) is greater than that inferred in a homogeneous model.Indeed, since the gravity contribution ∆gM is negative, the observed positive gravity change must be due to a larger density for the new mass.Also MOGI-5 (low rigidity and low density caldera) requires high density in spite of the positive ∆g M 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/m 3 .This value is much higher  II.Based on the evaluation of gravity residuals and the assumption of a pressurized point-source, Berrino et al. (1984) and Berrino (1994) suggested that a subsurface mass increase of 2 • 10 11 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 computed density is ∼2000 kg/m 3 while the mass input is 3-4 times the estimated value.On the other side, the mass calculated for the sill model is a factor 4 less than the reference value.

Conclusions
We develop a numerical technique to calculate 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 heterogeneities.Very few analytical solutions are available in literature (e.g., sphere, sill, spheroid in elastic media and point-sources in viscoelastic media), while the FE models may be characterized arbitrarily.The gravity calculations distinguish between the contribution of the medium and of the mass intrusion.Very different contributions come from the deformation of the medium, according to source shape.We separate the effect of the displacement of layers interfaces (∆gL) and of the medium compressibility (∆gV).We consider the gravity changes occurred at CF caldera during the 1982-1984 unrest 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 affecting the observations.It must be noted that the large deformation is well outside the elastic limit of any rock material.We show that nonspherical sources such as sill and spheroid yield positive or negative gravity changes (respectively), without input of new mass.The density values computed for these models are very different, changing by a factor 8 between extreme end models.The results clearly point out the importance of the source shape and of the characteristics of the medium in gravity calculation and source density estimation.
The proposed models, even if not yet oriented yet to data inversions, give hints in understanding the effects of considering the realistic structure of the studied area in gravity computations.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.Density 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 extracted from deep drilling seem the only plausible 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 procedure to retrieve the shape of a point-like deformation 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): including gravity data gives further constraints to restrict the range of acceptable models.

Fig. 1 .
Fig. 1.Map of Campi Flegrei caldera showing the outer rim (open triangles) and the inner caldera (full triangles), the leveling lines (dotted) were surveyed during 1982-1984 unrest episode, while the squares represents the gravity stations active in that period.
Fig. 2a-c.a) gravity variations and elevation change observed during 1981-2001 at Serapeo roman market (Pozzuoli), at the center of CF; b) gravity due to the ground uplift (free-air correction); c) residual gravity due to the difference between the observed and the free-air corrected gravity.

Fig
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 integrand 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).

Fig. 5 .
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.

Fig. 6 .
Fig. 6.Contributions and medium gravity change due to an inflating spherical pressure source in a homogeneous 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.

Fig
Fig. 7a,b.Gravity terms normalized to the maximum uplift uz max , and vertical displacement in a homogeneous medium due to (a) vertically elongated spheroidal source, (b) horizontal sill (in this case ∆gV 0 and ∆gM overlaps ∆gL).

Fig
Fig. 10a,b.Gravity terms normalized to the maximum uplift uz max , and vertical displacement in an inelastic medium: a) all plastic, b) plastic caldera and remaining elastic.

Fig. 11 .
Fig. 11.Density calculations versus volume variations of the proposed models.

Table II .
Volume variation and density calculated for each model (characteristics reported in tableI) due to input of new mass within the cavity.