Finite element modeling of ground deformation and gravity field at Mt . Etna

An elastic 3-D axi-symmetric model based on Finite Element Method (FEM) is proposed to compute ground deformation and gravity changes caused by overpressure sources in volcanic areas. The numerical computations are focused on the modeling of a complex description of Mt Etna in order to evaluate the effect of topography, medium heterogeneities and source geometries. Both ground deformation and gravity changes are investigated by solving a coupled numerical problem considering a simplified ground surface profile and a multi-layered crustal structure inferred from seismic tomography. The role of the source geometry is also explored taking into account spherical and ellipsoidal volumetric sources. The comparison between numerical results and those predicted by analytical solutions disclosed significant discrepancies. These differences constrain the applicability of simple spherical source and homogeneous half-space hypotheses, which are usually implicitly assumed when analytical solutions are applied. Mailing address: Dr. Gilda Currenti, Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Catania, Piazza Roma 2, 95125 Catania, Italy; e-mail: currenti@ct.ingv.it Vol51,1,2008_DelNegro 16-02-2009 21:28 Pagina 105


Introduction
Microgravity and geodetic observations have proven to be useful methods for monitoring volcanic activity and for the quantitative evaluation of the geophysical processes preceding and accompanying volcanic unrest (Carbone et al., 2003;Gottsmann et al., 2006).Ground deformation studies provide insight on volume changes in the magma reservoir and the dynamics of dike intrusion processes (Voight et al., 1998;Battaglia et al., 2003;Murase et al., 2006).However, deformation data alone are not able to properly constrain the mass of the intrusions.Geodetic studies need to be supported by gravity observations in order to infer the density of the intrusive body and better define the volcanic source (Carbone et al., 2006).
Attempts at modeling gravity changes expected to accompany crustal deformation often involve a great deal of effort due to the complexity of the problem.A series of analytical solutions for modeling ground deformation and gravity anomalies due to volcanic pressure sources have been devised and are widely used in literature (Okubo, 1992;Jousset et al., 2003).These models estimate the intrusive mass parameters (radius, depth, density, volume and pressure changes) from the geophysical variations observed at the ground surface.Most of the analytical formulations for modeling inflation and deflation episodes are based on the assumption of a pressurized magma chamber embedded in a homogeneous elastic half-space medium.These models describe the effects caused by sources with a specific geometry such as spheres (Hagiwara, 1977), ellipsoids (Battaglia and Segall, 2004) or rectangular prisms (Okubo and Watanabe, 1989).Analytical elastic models are quite attractive because of their straightforward formulation.However, volcanic areas are usually characterized by severe heterogeneities and consistent topography that are responsible for significant effects.Moreover, the simple assumption of spherical geometry for pressure sources can affect the estimate of depth and intrusion density, leading to a misinterpretation of the geophysical observations (Battaglia and Segall, 2004).Initial attempts to estimate the effect of heterogeneity were made solving a coupled self-gravitational model of magma intrusion in an elastic-gravitational layered medium (Rundle, 1980;Fernandez et al. 1999).These semi-analytical solutions assessed the effects of both elastic and density stratifications on gravity changes (Tiampo et al., 2004).Although numerical models over the past decade have emphasized that topography engenders perturbations on the surface deformation (Williams andWadge, 1998, 2000;Cayol and Cornet, 1998;Lungarini et al., 2005), few studies have been carried out to assess how topography together with source geometry can influence gravity anomalies.
To fill these gaps we propose a numerical procedure based on Finite Element Method (FEM) to evaluate both elevation and gravity changes expected in volcanic region.The numerical computations are focused on the modeling of a more complex description of Mt Etna.A 3D axi-symmetric model is developed assuming a multilayered medium for the subsurface structure and a conic shape for the volcano topography.Several tests are also carried out to appraise the influence of the volcano edifice on deformation and gravity fields.We also compare analytical and numerical solutions to estimate the inaccuracy that could be caused by neglecting medium heterogeneities and topography.Moreover the source geometry effect is investigated separately considering spherical and ellipsoidal sources in a homogeneous halfspace.Then, a comprehensive model is realized to understand how this more complex framework could affect the deformation and gravity field caused by likely overpressure source.

Numerical model
Generally, the gravity changes caused by pressure sources are made up by four different contributions (2.1) where δg0 represents the «free air» gravity change accompanying the uplift of the observation site; the δg1 is due to the displacement of density boundaries in heterogeneous media; δg2 denotes the effect of chamber wall expansion coming from pressurization effects which implies replacement of mass with gas or magma; and δg3 is the gravity change produced by density variations in the surrounding medium.
As for the first term in eq. ( 2.1), free-air correction only adjusts for the elevation of the observation point that, in first approximation, is given by where γ is the free-air gravity change rate (γ=308.6µGal/m) and δh the elevation change.
As for the other three terms in eq.(2.1), the gravity anomaly δgi, with i = 1, 2, and 3 indicating the different contributions, is related to the gravitational potential ϕg i as (2.3) where ϕg i can be inferred through the following Poisson's differential equations with appropriate boundary conditions (Cai and Wang, 2005). (2.4) These equations relate the gravity potential to the density variations caused by the subsurface mass redistribution: G denotes the universal gravitational constant, u is the displacement field, ρ 0 is the embedding medium density, ρmi is the source density before the magmatic intrusion and ρmf is the source density after the intrusion.
The term in eq.(2.4) is related to the density change due to the displacement of density boundaries in heterogeneous media.The δg2 gravity contribution (eq.(2.5)) depends on: i) the displacement of the source boundaries which implies replacement of surrounding mass (δg2∆V) and ii) the input of new mass inside the source volume (δg2V).Finally the δg3 term in eq.(2.6) takes into account the source rate of expansion arising from compressibility of the surrounding medium (Bonafede and Mazzanti, 1998).Therefore, it emerges that in order to evaluate these gravity contributions the computation of the displacement field at depth is required (Charco et al., 2006).Since gravity changes are strictly related to deformation of the elastic medium, the deformation field and the gravity anomalies produced by volcanic sources need to be modeled jointly.
We compare analytical and numerical results for all the gravity terms separately with the aim to quantify the inaccuracy introduced by the half-space assumption.In order to include the effect caused by topography, elastic and density medium heterogeneities and source geometry, we chose the finite element technique, suitable to solve both the elastostatic equation and Poisson's differential equations (eqs.(2.4) to (2.6)).Computations were carried out using the commercial software COMSOL Multiphysics, version 3.2.In the following we tackle this coupled numerical problem computing the displacement field and the gravity anomaly to estimate the effects of δg1, δg2 and δg3 terms on the total gravity variation.
We solve the model equations in two steps: i) the elastostatic equation in order to compute the elastic deformation field and its derivatives, ii) the coupled Poisson's gravity equation for gravity potential field to derive the gravity changes described above.

Benchmark test
The numerical model needs to set some parameters that could affect the accuracy of the solution.The size of the computational domain and the size of the finite elements are to be accurately chosen.The domain size is important because of the assignment of the boundary con-ditions.As for finite elements size, the meshing operation is a fundamental step: the smaller the elements, the more precise the solution.However, if too small elements are used, the number of nodes in which the equations are to be solved increases and computation becomes heavy.To properly set these numerical parameters, benchmark tests are carried out to calculate both crustal deformation and the gravity contributions.
A computational domain of a 100×100 km is considered for the deformation field calculations.As for boundary conditions, displacements are fixed to zero at the bottom and the lateral boundaries of the domain.These boundaries are located far away from the source so as not to measurably affect the results.The upper boundary is traction free and represents the ground surface.In order to solve Poisson's equation (eqs.(2.4) to (2.6)) the potential or potential derivative are to be assigned at the boundaries of the domain (Dirichlet or Neumann boundary conditions).The gravity potential is generally unknown on the external boundaries but it vanishes to zero at infinite distance (Zhang et al., 2004).The computational domain of the numerical problem is extended up to a distance of 100 km and null boundary conditions are imposed for the potential.In this way, the domain is big enough that the zero potential assumption does not affect the solution in the interested area.Neumann boundary conditions are instead imposed along the symmetry axis to satisfy the continuity of the gravity potential.
The deformation and the gravity changes are modeled for a spherical source, with an overpressure of 100 MPa and radius 0.5 km, buried at a depth of 5 km in a half-space homogeneous medium characterized by a density value ρ0 of 2500 kg/m 3 and a Young modulus of 62.5 GPa.Following this framework, the numerical solutions are compared with the analytical ones.The analytical solutions of the elastostatic and gravity problems are obtained using the simple and common Mogi model (Mogi, 1958;Hagiwara, 1977).The gravity changes δg1, δg2 and δg3 caused by the expansion of a spherical source embedded in a homogeneous Poisson's medium (λ=µ) are given by the fol-lowing analytical expressions (Hagiwara, 1977): (3.1) where ρmi (ρmf) and Vi (Vf) are the average source density and the volume source before (after) the inflation, ∆V=Vf −Vi is the volume change (Mogi, 1958) and r is the horizontal distance from the source.A good tradeoff between the computational domain and the computational load is found in a way to assure a good correspondence between analytical and numerical solutions (fig.1).It is also verified that whenever the source inflates without addition of new mass the overall gravity change δg1+δg2+δg3 vanishes identically (Walsh and Rice, 1979).The computations are carried out assuming a density change from ρmi=2400 kg/m 3 to ρmf=2800 kg/m 3 .Initially, we study the perturbations caused by the structural multilayering of the medium.Successively, we take into account simple topographic profiles to evaluate the influence of the topography on the solutions.Finally, we evaluate the effect due to the geometry of the pressurized volume.

Heterogeneity effect
A stratified medium for the subsurface structure of Mt Etna is considered assigning to each layer the elastic and density properties derived from seismic tomography data (Chiarabba et al., 2000).The density values are defined using empirical relationships between elastic wavespeeds and crust density (Birch, 1961;Brocher, 2005).As for the elastic parameters, if a Poisson medium (Poisson ratio equal to 0.25) is assumed, Young's modulus, E, is related to crust density, ρ0, and wavespeed, Vp, through the following relation: Starting from these assumptions, we consider a multi-layered crustal structure with six horizontal layers whose elastic and density parameters are reported in table I.The range of density values varies from 2000 kg/m 3 to 3500 kg/m 3 as the depth increases and they are in agreement with the crustal structure of Mt Etna inferred from compositional studies (Corsaro and Pompilio 2003).The elastic Young's modulus varies in the range from 11.5 GPa at shallower depth to 133 GPa at higher depth.
The gravity contributions δg1, δg2 and δg3 are evaluated for a pressurized spherical source with radius of 1 km and an overpressure of 100 MPa located at three different depths: (A) 3 km, (B) 6.5 km and (C) 9 km.The numerical solution of the multilayered model is compared with the analytical formulation in the case of a homogeneous half-space medium (ρ0=2500; E= =62.5 GPa).The elastic and density heterogeneities affect both the δh and the δg contributions depending on the source depth (fig.2).For the shallowest source (A), the δg1 ampli-   tude in the multilayered model at the center of the source is triple the homogeneous one.At larger horizontal distance from the source (>10 km), the heterogeneity does not give significant effects.For the source at medium depth (B) a less signal enhancement is obtained, while for the deepest source (C) quite similar results to the analytical solution are achieved.In the latter case, a discrepancy less than 1 µGal between the numerical and analytical solutions is obtained over the entire computational domain.
The δg2 term plays an important role in the discrimination between processes at constant mass and deformation mechanisms with mass intrusion.We perform numerical computations on both cases.When no mass is added inside the source, the law of mass conservation yields to where ρmf (ρmi) is the density and Vf (Vi) is the volume after (before) the expansion.For a spherical source, the volume change depends on the pressure change as where P is the source pressure, a is source radius and µ is shear modulus.Therefore, the density after the expansion can be inferred from (4.4) and it depends on the elastic medium parame- ters.After the expansion all three sources reach the same final volume and the same final density in the analytical model because of the homogeneous assumption.Instead, in the numerical model the volume expansion and the final density depend on the shear modulus of the layer in which the source is embedded (eq.(4.4)).In fig. 3 it is clear that the δg2V contribute diverges from the analytical one especially for the source (A), where the misfit is about 30 µGal, while for the other two sources there are no meaningful differences.The δg2∆V term shows a similar behavior, even if the discrepancies with respect to the analytical solution is less than 3 µGal for the sources (A).A reverse in sign is found for the sources (B) and (C) because after expansion the density contrast between the surrounding rock and the source becomes negative.Indeed, the amplitude of these gravity anomaly is too small to be taken into account, though this contribution could be remarkable whenever the sources would have larger dimensions.
When the source expansion is accompanied by mass intrusion, the δg2V contribution depends on the difference between the average density of the source before and after the inflation and on the distance between the observa-Fig.3. The δg2∆V and δg2V terms when no mass is added inside the source using the heterogeneity model (dashed lines) for a source at three different depths: 3 km (blue), 6.5 km (green), 9 km (red).The analytical solutions for homogeneous medium are also shown (solid lines).tion point and the source.Therefore, this contribution is not affected by the medium heterogeneities.As a result, in the following we numerically compute only the δg2∆V term since no differences were found between the numerical and the analytical solution for the δg2V term (eq.(3.2)).A density of 2800 kg/m 3 after the inflation is assumed.Results show that the maximum discrepancy in the δg2∆V term is once again for the shallower source (fig.2).
As for the δg3 contribution, the pressurized source (A) in the multilayered medium gives rise to a negative gravity anomaly of about −35 µGal above the center of the source.Similar sources in a homogeneous medium bring only about −11 µGal.Figure 2 shows significant differences between the heterogeneous and homogeneous model mainly on the ground surface above the source.The (B) and (C) heterogeneous models generate gravity anomalies quite similar to the analytical ones (fig.2).Indeed, the misfit with respect to the homogeneous model is well within the measurement error for gravity data acquired in volcanic areas.Therefore, for all the gravity contributions, significant deviations from the homogeneous case are obtained especially for source (A).We can conclude that the shallower the source is the higher the discrepancy between the multilayered model and the homogeneous one.When the medium is elastic, the deformation pattern and the gravity changes δg1, δg2∆V and δg3 depend almost entirely on the average rigidity and the density of the medium surrounding the source.In fact, the shallower source is surrounded by a lower rigidity medium that leads to a higher displacement around the pressure source.The result is an increase in the density variation because of the enhancement of the u displacement and the ∇⋅u rate of expansion (eqs.(2.4) to (2.6)) in the region surrounding the source.

Topography effect
In order to quantify the effect of topography on gravity anomalies, we perform several tests by adding a topographic relief to the previous multilayered models.We model Mt.Etna topography using the a simplified ground surface (fig.4).The volcano edifice is represented as a cone 3300 m high and 20 km in radius with an average slope of the flanks of about 10°.As described above, the computations are performed using the pressure source at three different depths.To estimate the effects of topography on gravity anomalies a comparison with the analytical solution is carried out.In the analitic model we choose a reference surface at 1500 m that represents the average altitude of the topographic relief.For all three pressure sources, the δg1 contribution is enhanced with respect to the analytical solutions.What is more important is the change in shape observed in the proximity of the volcano summit, where a sharper decrease is obtained (fig.5).At a horizontal distance of about 10 km from the summit the numerical results are quite similar to the analytical ones.
Following the same procedure as for the multilayered model, the δg2 contribution is computed in the case of mass intrusion and constant mass.When no mass intrudes, the δg2V and δg2∆V gravity terms are small because the effective distance between the source and the surface is enhanced.The perturbations intro- duced by the numerical model have no importance especially in the δg2∆V term (fig.6).For a Mogi model, when the source inflates with mass intrusion, the topography effect on the δg2V term can be simply accounted for replacing the source depth z by zl, where zl= z + f with f the elevation above the reference level (Charco et al. 2007).So even for the topography effect, we calculate only the δg2∆V term.The δg2∆V term shows variations less than 2 µGal and therefore is negligible for all three sources.
In addition the δg3 gravity term does not show significant shape deviations from the analytical results when topography and heterogeneity features are simultaneously accounted for.The enhancement in the δg3 contribution depends on the depth of the source.For the shallower source this contribution is almost doubled, while for the deeper source the increase is less evident.It is worth noting that, summing up all the terms δg1+δg2+δg3 in the case of constant mass process, the overall gravity change does not vanish because of medium heterogeneities and topography.Our findings, shown in figs. 2 and 5, show that the analytical models underestimate elevation changes and gravity variations.Therefore, if analytical models are used to estimate the volcanic source parameters on the basis of microgravity and geodetic data observed at the ground surface, the Fig. 5. Numerical (dashed lines) and homogeneous half-space solutions (solid lines) for the elevation changes and the δg1, δg2∆V and δg3 gravity contributions, including the topography effect.The increase in the gravity field depends on the depth of the source: 3 km (blue), 6.5 km (green), 9 km (red).
pressure change would be overestimated, while the source depth would be underestimated.

Source geometry effect
Most of the formulations for modeling inflation and deflation episodes assume simple geometry sources, commonly a point source, which well approximates a spherical pressurized cavity when the source depth is greater than twice the cavity diameter (McTigue, 1987).This simple source approximation has been adequate to interpret many geodetic and gravity data in volcanic areas (Dvorak and Dzurisin, 1997;Fernandez et al., 1999).The model geometry used to invert the data might bias the depth and volume estimate of the inflation source.In fact, it can be seen that spherical model systematically yields a deeper source than ellipsoidal one and therefore requires a larger mass to obtain the same gravity changes (Battaglia et al., 2003).Hence, the type of the assumed source can lead to an inaccurate estimate of the depth and thus density of the intrusion.In order to investigate the role of the source geometry, we compare the vertical displacement and the gravity field produced by three different sources: a spherical source (SP), a vertical prolate spheroid source (VPS, dykelike shape) and a horizontal prolate spheroid source (HPS, sill-like shape).For all three sources the geometric parameters are set to ensure the same volume.We fix the radius for the spherical cavity equal to 500 m.Since the two equatorial radii (a, b) of the spheroids are equal, the choice is in term of the ratio between the equatorial and the polar (c) radii that is fixed as four times the equatorial radii for both the spheroids (a=b=314.98m, c= 1259.92m for VPS, a=b=1259.92m, c=314.98 m for HPS).Previous results have showed that major effects are observed for the shallower sources; with this in mind the numerical computations are performed positioning the source at a depth of 1.5 km.
Firstly, we numerically model the vertical deformation and the δg1, δg2 and δg3 gravity terms in a Poissonian homogeneous half-space (Young's modulus =62.5 GPa and density =2500 kg/m 3 ) to separately appraise the geometry effect.Indeed, the spherical source and the Fig. 6.Numerical (dashed lines) and homogeneous half-space solutions (solid lines) for the δg2∆V and δg2V gravity contributions whether no mass intrudes into the source, including the topography effect.The increase in the gravity field depends on the depth of the source: 3 km (blue), 6.5 km (green), 9 km (red).vertical ellipsoid behave in a similar way for both deformation and gravity terms (fig.7).The HPS source causes a vertical displacement, δh, of about 2 m against the 15 cm causes by the SP and VPS sources.It is worth noting that the δg1 term for the HPS source differs by 200 µGal respect to the VPS and SP sources.These differences are considerable even at a horizontal distance of 2-3 km from the source.Comparable results are obtained for the δg3 gravity term although the discrepancy is reduced to about 12 µGal.The VPS model effect decays faster than the HPS one in this case.The δg2V term is not reported since it coincides with the analytical solutions for ellipsoidal sources devised by Clark et al. (1986).Furthermore, the geometry effect for the δg2∆V term could be neglected, because the expected gravity changes are within a range of 2 µGal for all the sources considered (fig.7).
Summing up, the HPS yields gravity changes that differ by about an order of magnitude with respect to those generated by the other two sources even if they are subjected to equal volume change and the top of the HPS Fig. 7. Elevation and gravity changes for three source shapes: the SP source (solid line), the VPS source (dashed line) and the HPS source (dotten line).The numerical model assumes an homogeneous elastic half space.The geometric parameters were chosen for all the sources in a way to assure the same volume.The source depth was fixed to 1.5 km.source is lower.This is caused by the different force line distributions acting on the spheroid surface in the HPS and the VPS models.Indeed, the line forces in the HPS are for the most part orthogonal to the ground surface, while in the VPS they are mainly parallel: the greater the z-component of the surface forcesis, the greater the intensity of gravity terms.Similar results were obtained for ground deformation by Russo and Giberti (2004).To better understand if the relations between δh and the δg1 and δg3 gravity changes are affected by the source geometry, the gravity terms are referred to the vertical displacement for all three sources.The δg1/δh ratio remains constant independently from the geometry, whereas the δg3/δh ratio reveals a particular feature.The two prolate spheroid sources exhibit a non linear pattern that differs from the well-known linear relationship in the SP model (eq.(3.3)).In proximity of the volcano summit the ratio between δg3 and δh is strongly dependent on the source geometry (fig.8).A small variation in δh leads changes in the δg3 gravity contribution that are higher for the VPS source than for the HPS one.
Finally, we analyze the all-in effect in the displacement and gravity field when source geometry, medium heterogeneities and topography are simultaneously accounted for.Three   sources, a spherical one and two prolate spheroids (a=b=396.8m, c= 793.7 m), are located at 0 km depth in the above described multilayered axi-symmetric model with a topographic profile (fig.9).Topography dampens the source geometry effect because the effective distance to the source is increased (fig.10).As before, the HPS gravity anomalies strongly differ from VPS and SP sources, by about 1 m in the vertical displacement, 50 µGal in the δg1, 2 µGal in the δg2∆V and 15 µGal in the δg3 terms.Moreover, the δg1 term shows a change in shape with a reverse in sign at about 3 km from the summit.This could be generated from the gravitational attraction of the gain of density increase above the observation points due to the elevation changes of the local topography.In fact, this negative anomaly is observed in the lower altitude points located near the volcano summit.
Because of the topography, the δg2V needs to be numerically evaluated.Assuming a density change from ρmi=2400 kg/m 3 to ρmf=2800 kg/m 3 , the δg2V gravity term exhibits a discrepancy of about 150 µGal when the geometry changes from SP to HPS and about 100 µGal from VPS to SP.

Discussion and conclusions
The finite element model is applied to evaluate deformation and gravity changes produced by pressurized volcanic sources.The FEM method allows us to investigate three factors that can help to obtain a more realistic picture of the volcanic framework and of the intrusive body: topography, elastic heterogeneities and source geometries.The results presented are obtained developing a 3D axi-symmetric numerical model, in which heterogeneities of country rocks and the shape and location of the overpressure source play an essential role.Among the four gravity contributions generated by pressure sources, we do not report on the gravity term due to free air effect (δg0), since gravity changes are usually corrected, before the interpretation, for the effect of uplift (the free-air correction) using leveling data.We take into account not only the elastic heterogeneities of the medium but also the density distribution of subsurface structures.Therefore, the δg1, which usually accounts for only the Bouguer anomaly, is computed considering also the displacements of the subsurface density layers.Numerical solutions for the δg1, δg2 and δg3 terms, computed by the finite element method, have been compared to analytical results for a homogeneous half-space medium.To appraise the influence of source geometry on the solution, we have evaluated elevation and gravity changes for spherical, vertical prolate spheroid (dyke-like shape) and horizontal prolate spheroid (sill-like shape) sources.
Heterogeneity, topography and source geometry engender deviations from analytical results in the deformation and gravity changes produced by pressurized sources under elastic conditions.However, such perturbations are more evident in the presence of severe heterogeneities and steeper topography, i.e. in the volcano summit.The elastic heterogeneity only affects the amplitudes of anomalies, while the topography strongly alters both the magnitude and the shape especially in the δg1 contribution.In particular, our results evidenced that the analytical model underestimates elevation changes and gravity variations.Moreover a critical point is the choice of the source geometry used to interpret geodetic and gravity data.Indeed the elevation and gravity changes for the spherical and ellipsoidal sources exhibit significant differences that are strongly emphasized because of the different force line distributions acting on cavity surfaces.Therefore, the simple spherical geometry assumption can lead to a wrong estimate of source parameters (depth, mass and density of the intrusion) in inverse problem.
Our findings highlighted two main points.Firstly, changes in the gravity field cannot be interpreted only in term of additional mass input disregarding the ground deformation of the surrounding rocks.Hence, in order to obtain a reliable estimate of the depth and density of the intrusion, geodetic and gravity data, which independently reflect the state of volcano, need to be jointly modeled.Secondly, standard analytical models that neglect the complexities associated with morphology and medium properties of volcanic edifice, could provide an inaccurate estimate of the ground deformation and gravity anomalies expected at the ground surface.This can lead to a misinterpretation of the geophysical observations carried out in volcanic areas.The application of finite element methods, instead, allows for a more accurate modeling procedure, which might provide sensible insights into volcanic source definition.

Fig. 1 .
Fig. 1.Comparison between the analytical solutions (solid line) and the finite element results (circles) for the δg1 (blue), δg2∆V (red) and δg3 (green) for a spherical source.The overall gravity change vanishes identically (black circle) as demonstrated analytically (black solid line) by Walsh and Rice (1979).

Fig. 2 .
Fig.2.Elevation and gravity changes using the heterogeneity model (dashed lines) for a source at three different depths: 3 km (blue), 6.5 km (green), 9 km (red).The analytical solutions for homogeneous medium are also shown (solid lines).

Fig. 4 .
Fig. 4. Model geometry used to estimate the influence of topography on surface deformation and gravity field.The features in the picture are not in scale.

Fig. 8 .
Fig.8.Ratios between δg3 gravity terms and δh elevation changes using the numerical model for three source shapes: the SP source (solid line), the VPS source (dashed line) and the HPS source (dotted line).

Fig. 9 .
Fig. 9. Model geometry used to estimate the influence of topography and hetereogeneities on surface deformation and gravity field for the three geometry sources: HPS, VPS and SP.The arrows represent force line distributions acting on the source surface.The features in the picture are not in scale.

Fig. 10 .
Fig. 10.Numerical solutions for the elevation changes and the δg1, δg2∆V and δg3 gravity contributions for three source shapes: the SP source (solid line), the VPS source (dashed line) and the HPS source (dotted line).The numerical model takes into account elastic and density heterogeneities of the medium and topography.

Table I .
Medium properties of the multilayered model.