# A strategy to explore the topography-driven distortions in the tilt field induced by a spherical pressure source: the case of Mt Etna

François Beauducel, Daniele Carbone

► To cite this version:

François Beauducel, Daniele Carbone. A strategy to explore the topography-driven distortions in the tilt field induced by a spherical pressure source: the case of Mt Etna. Geophysical Journal International, 2015, 201, pp.1471-1481. 10.1093/gji/ggv076. insu-03579969

HAL Id: insu-03579969

https://insu.hal.science/insu-03579969

Submitted on 18 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A strategy to explore the topography-driven distortions in the tilt field induced by a spherical pressure source: the case of Mt Etna

François Beauducel<fnref n="1" /> and Daniele Carbone<fnref n="2" />

Accepted 2015 February 16. Received 2015 February 13; in original form 2014 July 8

## SUMMARY

We present a strategy to thoroughly investigate the effects of prominent topography on the surface tilt due to a spherical pressure source. We use Etna's topography as a case of study and, for different source positions, we compare the tilt fields calculated through (i) a 3-D boundary element method and (ii) analytical half-space solutions. We systematically determine (i) the source positions leading to the strongest tilt misfits when numerical and analytical results are compared and (ii) the surface areas where the strongest distortions in the tilt field are most likely to be observed. We also demonstrate that, under critical circumstances, in terms of respective positions of pressure source and observation points, results of inversion procedures aimed at retrieving the source parameters can be misleading, if tilt data are analysed using models that do not account for topography.

**Key words:** Numerical approximations and analysis; Transient deformation; Volcano monitoring.

## 1 INTRODUCTION

Measurements of surface deformation at active volcanoes can be used to study the dynamics of subsurface magma storages and to constrain the mechanisms leading to eruptions (Dvorak & Dzurisin 1997; Dzurisin 2003; Segall 2010). Volcanic deformation is often modelled using analytical solutions that offer a closed-form description of the source. Simplifications, such as the assumption that the crust is a homogenous, isotropic, elastic, semi-infinite medium with flat upper surface (elastic half-space), make analytical models tractable and easy to implement, but can lead to misleading volcanic interpretations (Lisowski 2006; Battaglia et al. 2008).

Anderson (1936) proposed analytical solutions to an isotropic point source of pressure, that can be used to approximate a spherical magma chamber (Mogi 1958). The analytical solution of Okada (1985) permits to calculate surface deformation from finite rectangular dislocations and can be used to model the effect of fault dislocations or dyke-like intrusions. In spite of the steep topography of many volcanic areas, these half-space solutions are still widely used (e.g. Aloisi et al. 2006; Anderson et al. 2010), especially because they are easy to use as forward models for inverse methods. Early attempts to account for topography-induced distortions of the observed deformation field were carried out using 2-D analytical solutions (Ishii & Takagi 1967; Harrison 1976; Meertens & Wahr 1986; McTigue & Segall 1988). Through comparison with results from a 3-D finite element model, Williams & Wadge (1998, 2000) explored the capabilities of semi-analytic solutions to account for the effect of topography on volcano deformation.

Topography-induced distortions in surface deformation are not only due to the different elevations of the observation points and the slope gradient between them. Rather, these distortions are due to the real 3-D topographic features of the volcanic edifice. This is particularly the case for tilt which is more sensitive to local effects. The use of boundary or finite elements methods allows to accurately and exhaustively account for the effect of topography on the surface deformation induced by a buried pressure source (e.g. Cayol & Cornet 1998; Beauducel & Cornet 1999; Fukushima et al. 2005; Currenti et al. 2007; Meo et al. 2008; Trasatti et al. 2008; Lyons et al. 2012; Masterlark et al. 2012; Charco & del Sastre 2014). However, complex modelling schemes cannot be always used, especially for real time volcano monitoring, where computation time is a critical issue. On the other hand, neglecting the topography may lead to erroneous conclusions on the characteristics of the deformation source (e.g. position and volume/pressure variation) and even to inversions of the sign of the pressure change at the source (Cayol & Cornet 1998; Beauducel & Cornet 1999). It is thus useful to a priori identify the circumstances where the most significant distortions in the deformation field can arise, if topographic effects are disregarded.

In the following, we describe a strategy to explore the topography-driven distortions in the tilt field induced by a spherical pressure source. We focus on tilt changes (and we disregard, for example vertical/horizontal displacements) because of their ability to image processes that develop over short-timescales (from minutes to hours; Ohminato et al. 1998; Voight et al. 1998; Anderson et al. 2010; Montgomery-Brown et al. 2010; Peltier et al. 2011) and that are of

<footnote n="1">¹Institut de physique du globe de Paris, Sorbonne Paris Cité, Univ. Paris Diderot, UMR 7154 CNRS, France. E-mail: beauducel@ipgp.fr</footnote>

<footnote n="2">²Instituto Nazionale di Geofisica e Vulcanologia, Sezione di Catania - Osservatorio Etneo, Italy</footnote>
primarily interest for real-time monitoring. Furthermore, short-base
tiltmeters are expected to be more sensitive to topographic effects
than, for example, GPS instruments, since they measure the spatial
derivative of the displacement field (Harrison 1976; van Driel et al.
2012).

Our exploration strategy is implemented by comparing the re-
sults from a 3-D numerical method with results from analytical
solutions for a homogeneous half-space. Eventually, we determine
(i) the source positions at which the maximum differences between
the results from numerical and analytical models occur and (ii) the
surface areas where the strongest misfits are more likely to occur.
We also show and discuss how the distortions arising from Etna's
topography on the tilt field induced by a pressure source may af-
fect the inversion of the source parameters. The choice of Etna's
topography as a case of study is motivated by the strong asym-
metry of its volcanic edifice that results from a complex eruptive
history and an active tectonic environment (Froger et al. 2001). The
seaward-facing eastern flank is deeply dissected by normal faults
and presents significant signs of instability. The most striking of
Etna's features is the Valle del Bove (hereafter VdB; Fig. 1a), a
deep horse-shoe shaped caldera resulting from repetitive collapse
sequences in the last 100–150 kyr (Calvari et al. 2004).

## 2 METHODOLOGY

To investigate the effect of Etna's topography on the tilt field in-
duced by a spherical pressure source, we compare the corresponding
tilt fields calculated through the mixed boundary element method
(MBEM) developed by Cayol & Cornet (1997), the analytical so-
lution (MOGI) developed by Mogi (1958) and the 'depth-varying'
analytical solution (DVAS) proposed by Williams & Wadge (1998,
2000).

The MBEM is a 3-D scheme that allows the analysis of elastic
deformation fields, taking into account topography, medium discon-
tinuities and pressure sources of any shape. 3-D boundary element
methods (Cianetti et al. 2012; Privitera et al. 2012) do not require
complex meshing or long computation time. MBEM is a combi-
nation of two boundary element methods: (i) the direct method,
based on Betti's reciprocal theorem (Rizzo 1967) and the solution
of Kevin's problem of a point force in an infinite body, and (ii)
the displacement discontinuity method (Crouch 1976). The direct
method allows accurate and fast modelling of structures such as
topography, cavities or reservoirs, whereas the displacement dis-
continuity method is suitable for fractures. The edifice is assumed
to be linearly elastic, homogeneous and isotropic. The prescribed
boundary conditions are tractions; they represent perturbations of
an initial state of stress, and are null on the ground and equal to
overpressure on deformation sources. The precision of MBEM was
carefully tested by Cayol & Cornet (1997) and the method has since
been successfully applied to several studies (e.g. Cayol & Cornet
1998; Beauducel et al. 2000; Cayol et al. 2000; Dieterich et al.
2000; Fukushima et al. 2005).

Here, Etna's topography, derived from a DEM with 10 m resolu-
tion, is described by a mesh composed of about 10 000 Delaunay
triangles. The size of the mesh elements increases with distance
from the centre of the mesh (summit zone of the volcano, min-
imum size of the elements = 0.2 km) to the periphery (size of
the elements = 2–3 km at the far edges of the mesh). Indeed, as
explained by Cayol & Cornet (1997), MBEM's results are more
accurate when the surface topography is discretized onto a graded
mesh (finer elements in the most perturbed area) than onto a regular

mesh. The graded mesh also requires fewer elements and is thus
computationally less expensive. The mesh we utilize here has an
area of 40 × 40 km and it is large enough to avoid border effects
(Fig. 1b). The source is modeled by a 100 m-diameter spherical
cavity represented by 1000 elements with constant pressure.

MOGI is based on the analytical solution of ground deforma-
tion induced by an isotropic point source of pressure in an elastic
half-space (Mogi 1958). MOGI solution is valid when the size of
the magma reservoir is very small compared to its depth. DVAS
differs from MOGI in being a quasi-analytic method where the ver-
tical distance from the source to the free surface varies to account
for the elevation of each computation point. This makes a simple
first-order account for the topography effect (Chaussard & Amelung
2012). Comparison between the results of MBEM and MOGI re-
veals the distortion in the tilt field arising when the real topography
is approximated by a flat surface at a constant elevation. Conversely,
comparison between the results of MBEM and DVAS reveals the
distortions in the tilt field arising only from the local shape of the
surface, since changes in the vertical distance from the source to
the free surface are accounted for by DVAS. To perform exhaustive
comparison between the results of the three methods, we use a trial
source inducing a constant pressure variation. The trial source is
placed at the (x, y, z) nodes of a 10 × 10 × 7 km grid, (hereafter re-
ferred to as 'source grid'; Fig. 2) with a fixed node-spacing of 1 km
(green dots in Fig. 1b). The source grid is centred on the location
of the summit craters (500 km east and 4178 km north; in UTM-
WGS84 coordinates), with top and bottom at elevations of 1 km a.s.l.
and 6 km b.s.l., respectively. For each source position, the tilt field at
the surface is computed using MBEM, MOGI (reference surface at
an elevation of 1.6 km a.s.l.) and DVAS. MBEM results are obtained
at the elements of the fine-mesh region (18.5 × 18.5 km area inside
the solid green square in Fig. 1b; size of the mesh elements equal to
0.2 and 0.4 km, in the inner and outer parts, respectively), by spatial
derivation of the displacement field. Results are then interpolated
at the (i, j) nodes of an uniformly spaced (100 m) grid covering the
same surface. In the following, we will refer this grid of 186 × 186
nodes as 'computational grid' (Fig. 2). MOGI and DVAS results are
obtained directly at the nodes of the computational grid.

It is worth pointing out that, since we consider only the misfit
between the results of analytical and numerical models due to the
effect of topography, we disregard possible medium heterogeneities
and we use the same elastic parameters (Young's modulus and Pois-
son's ratio) in all computations.

## 3 BENCHMARK OF MBEM AGAINST AN
EXACT ANALYTICAL SOLUTION

We benchmark the results predicted by MBEM through comparison
with Mogi's exact analytical solution. Since the latter is based on
the half space approximation, to perform the above comparison, we
assume that the top surface of the MBEM domain is flat. In practice,
we utilize the same mesh described above (Fig. 1b), with the only
difference that all the elements of the mesh are at the same elevation
(corresponding to the elevation of the reference surface of MOGI).

Fig. 3 shows corresponding radial and tangential tilt calculated at
all the nodes of the computational grid through the analytical (dotted
curve) and numerical (grey dots) schemes. Data are presented in
cylindrical coordinates, with the axis through the centre of the grid.
Tilt is induced by the spherical pressure source described in the
previous section (diameter = 100 m) located at 6 (Fig. 3a) and 1 km
(Fig. 3b) below the centre of the free surface. As for the deeper
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Topography of Mt Etna. The top-right inset shows the front of the Apennine-Maghrebian thrust belt and the Malta Escarpment fault zone, offshore eastern Sicily. (b) Sketch map showing the boundary element mesh of Mt Etna, that consists of about 10 000 Delaunay triangles. Green dots indicate the surface projection of the source grid. The solid green square encloses the fine-mesh region and the 100 m-spaced computational grid (see text for details).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2. 3-D Sketch of the source and computational grids.</figcaption></figure>

source, comparison between analytical and numerical predictions reveals a standard deviation of 1.5 and 0.7 per cent, for radial and tangential components, respectively. In the case of the shallower source, the standard deviation for radial and tangential components is 0.2 and 0.1 per cent, respectively. These small discrepancies are mainly attributed to numerical effects due to mesh discretization. Overall, the very good agreement between analytical and numerical predictions supports the accuracy of the MBEM results.

## 4 DESCRIPTION OF THE RESULTS

The strategy presented in the following subsections is aimed at describing (i) to which critical source positions in the source grid (Section 2; Figs 1b and 2) the strongest distortions in the surface tilt field, due to topographic effects, correspond and (ii) where on the surface the strongest distortions in the tilt field, due to topographic effects, are more likely to occur. The results of this twofold strategy thus permit to identify those circumstances, in terms of respective

positions of pressure source and observation points, where the most significant distortions in the tilt field induced by a spherical pressure source can arise, if topographic effects are disregarded.

As a common preliminary step of our twofold strategy (see the following subsections), we calculate, at each $(i, j)$ node of the computational grid (Fig. 1b) and for each $(x, y, z)$ position of the test source in the source grid (Fig. 2), the differences $m_{i,j}^{x,y,z}$ between (i) the amplitudes of the tilt vectors predicted through MBEM and MOGI models, (ii) the amplitudes of the tilt vectors predicted through MBEM and DVAS models, (iii) the directions of the tilt vectors predicted through numerical and analytical models. Note that, since the direction of the tilt vector does not depend on the vertical distance between source and observation point, the MBEM/MOGI and MBEM/DVAS differences in tilt direction, at the same observation point and for the same source positions, are identical. Eventually, a data set of about $9 \times 10^7$ differences between the predictions of numerical and analytical models is obtained.

### 4.1 Effect of source position

In the first step of our analysis we compute, for each given $(x, y, z)$ position of the test source in the source grid (Figs 1b and 2), the rms of the differences $m_{i,j}^{x,y,z}$ obtained at all the $(i, j)$ nodes of the computational grid (Fig. 2):

$$M_1^{x,y,z} = \sqrt{\frac{\sum_i \sum_j (m_{i,j}^{x,y,z})^2}{186^2}} \quad (1)$$

In general, the integrated $M_1$ misfits allow to identify the critical source positions to which strong distortions in the tilt field, due to topographic effects, are most likely to correspond.

From the above three categories of $m_{i,j}^{x,y,z}$ differences, we obtain three types of integrated $M_1$ misfits, namely: the MBEM/MOGI integrated amplitude misfit ($M_1M/M_{IAM}$); the MBEM/DVAS

<figure><img src="image_5.png" /><figcaption>Figure 3. Corresponding radial and tangential tilt calculated through analytical (MOGI; dotted curve) and numerical (MBEM; grey dots) schemes. To give a complete picture of the discrepancy between analytical and numerical predictions, data from all the nodes of the computational grid are presented in cylindrical coordinates, with the axis through the centre of the grid. The top surface of the MBEM domain is assumed to be flat and at the same elevation as that of the reference surface of MOGI. Tilt is induced by a spherical pressure source located at 6 (a) and 1 km (b) below the centre of the free surface. Data are normalized to the maximum amplitude of the radial component predicted by MOGI.</figcaption></figure>
integrated amplitude misfit ($M_1M/D\_IAM$); the numerical versus analytical integrated direction misfit ($M_1IDM$).

We show the integrated $M_1$ misfits (Fig. 4) in plan views corresponding to test sources at 6, 3, 1 and 0 km b.s.l. We choose these source depths as a fair compromise between the needs of adequately convey the relevant information, while avoiding redundancy. Results relative to source positions above the sea level are not shown since these positions are very close to the ground surface, especially in correspondence of the easternmost margin of the computational grid, where they would induce very strong tilt changes. Given the strong amplitude of these changes, even relatively small misfits between the predictions of analytical and numerical models would results in overwhelming effects and thus the tilt misfit maps would be dominated by local effects in correspondence of steep topographic features at lower elevations.

The maps of Fig. 4, are restricted to the area covered by the vertical projection of the source grid, while the integrated $M_1$ misfit for each source position is calculated over the whole computational grid. It is worth remarking that the maps in Fig. 4 do not represent tilt fields; rather, they represent integrated misfits for sets of source positions at different depths.

The values of the maximum misfits reported in Fig. 4 and Table 1 have no physical meaning as they result from integrating the differences between the effects predicted using either analytical or numerical models. Nevertheless, these values can be used to make comparisons between integrated $M_1$ misfits for pressure sources at different depths. The $M_1M/D\_IAM$s (left column in Fig. 4) mimic the surface elevation of the volcano and the resemblance becomes more obvious as the test source is placed at increasingly shallower depths. If one considers that the MOGI model is based on the

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Maps of integrated tilt misfits ($M_1$; see eq. 1) for source positions at different depths. Left column: MBEM/MOGI integrated amplitude misfit ($M_1M/D\_IAM$); middle column: MBEM/DVAS integrated amplitude misfit ($M_1M/D\_IAM$); right column: numerical/analytical integrated direction misfit ($M_1\_IDM$). Numbers reported on each panel represent relative maximum values of the colour scale (the minimum value is always zero; see Table 1 and the text for details).</figcaption></figure>
<table><caption>**Table 1.** Maximum misfit values $M_1$ and $M_2$ for different source depths, using a constant pressure source. Amplitude misfits are normalized to MBEM average amplitude, while direction misfits are expressed in degrees, relative to MBEM direction. Values correspond to the maximum of colour scale used in Figs 4 and 5.</caption><thead><tr><th rowspan="2">Source depth (km b.s.l.)</th><th colspan="3">Maximum misfit M<sub>1</sub></th><th colspan="3">Maximum misfit M<sub>2</sub></th></tr><tr><th>Ampl. MOGI (per cent)</th><th>Ampl. DVAS (per cent)</th><th>Direction (°)</th><th>Ampl. MOGI (per cent)</th><th>Ampl. DVAS (per cent)</th><th>Direction (°)</th></tr></thead><tbody><tr><td>0</td><td>2.6</td><td>0.38</td><td>4.4</td><td>73</td><td>17</td><td>2.6</td></tr><tr><td>1</td><td>1.1</td><td>0.21</td><td>5.7</td><td>31</td><td>9.2</td><td>2.2</td></tr><tr><td>2</td><td>0.64</td><td>0.18</td><td>6.9</td><td>19</td><td>6.6</td><td>2.6</td></tr><tr><td>3</td><td>0.44</td><td>0.16</td><td>8.6</td><td>14</td><td>5.7</td><td>3.0</td></tr><tr><td>4</td><td>0.33</td><td>0.14</td><td>10</td><td>11</td><td>5.4</td><td>3.4</td></tr><tr><td>5</td><td>0.27</td><td>0.13</td><td>12</td><td>9.6</td><td>5.4</td><td>4.0</td></tr><tr><td>6</td><td>0.23</td><td>0.13</td><td>14</td><td>8.2</td><td>5.4</td><td>4.3</td></tr></tbody></table>

assumption of a flat free surface at constant elevation, this result indicates that the effect induced by the different vertical distance between source and surface dominates over the effect induced by the local slope of the surface. Indeed, the $M_1M/D\_IAMs$ do not show the same correlation with surface elevation and are smaller (by a factor of between 2 and 7; see Table 1) than the corresponding $M_1M/D\_IAM$ (middle column in Fig. 4). In particular, for sources at greater depths, the strongest $M_1M/D\_IAMs$ are induced by sources placed below the summit of the volcano and below the VdB. At shallower depths the $M_1M/D\_IAMs$ due to sources located right below highly sloping features of the topography (e.g. margins of the VdB, NW upper flank of the volcano) become more important and, for sources placed at the shallowest depth of the chosen range (sea level), the effect of steep topographic features at lower elevations, thus closer to the test sources (e.g. southern margin of the VdB), become overwhelming.

The patterns of source positions corresponding to the strongest $M_1IDMs$ (maps in the right column of Fig. 4) mimic the topography of Etna, especially when the test sources are towards the deepest end of the chosen depth range. Stronger $M_1IDMs$ occur when the source is placed at deeper levels (Table 1). At shallow depths, the $M_1IDMs$ due to sources located right below highly sloping features of the topography become more important.

## 4.2 Effect of the position of the observation point

In the second step of our analysis, we calculate, at each single $(i, j)$ node of the computational grid, the rms of the differences $m_{i,j}^{x,y,z}$ obtained for all the source positions in a given horizontal layer of the source grid:

$$M_2^{i,j,z} = \sqrt{\frac{\sum_x \sum_y (m_{i,j}^{x,y,z})^2}{11^2}} \quad (2)$$

In general, the integrated $M_2$ misfits can be used to identify those areas where the strongest distortions in the tilt field, due to topographic effects, occur more frequently.

From the three categories of $m_{i,j}^{x,y,z}$ differences described at the beginning of Section 4, we obtain three types of integrated $M_2$ misfits, namely: the MBEM/MOGI integrated amplitude misfit ($M_2M/D\_IAM$); the MBEM/DVAS integrated amplitude misfit ($M_2M/D\_IAM$); the numerical versus analytical integrated direction misfit ($M_2IDM$).

As in Section 4.1, the integrated $M_2$ misfits are shown in plan views from 6, 3, 1 and 0 km b.s.l. (Fig. 5). As noted in the previous section, while the values of the maximum misfits reported in Fig. 5 and Table 1 have no physical meaning, they can nevertheless be used to make comparisons between integrated

$M_2$ misfits for pressure sources at different depths. The strongest $M_2M/D\_IAMs$ (left column in Fig. 5) occur more frequently in areas of higher elevation, that is where the largest differences between real and reference elevations arise. The normalized maximum misfit becomes systematically stronger with decreasing source depth (Table 1).

The strongest $M_2M/D\_IAMs$ (middle column in Fig. 5) occur more frequently in areas of steep topography (e.g. margins of the VdB and NW upper flank of the volcano). For sets of source positions at shallower depths, the strongest $M_2M/D\_IAMs$ misfits tend to concentrate around the highly sloping features of the topography at lower elevations. Also in this case, the normalized maximum misfit becomes stronger with decreasing source depth (Table 1). The stronger $M_2\_IDMs$ (right column in Fig. 5) are found more frequently on the vertical of steep topographic features. It is worth stressing that, especially for sources at shallower depths, there appear local spikes above the source positions in the $M_2\_IDM$ maps (see Fig. S1 in the Supporting Information). Indeed, since models with axisymmetric deformation about the vertical axis (e.g. MOGI and DVAS) predict opposite tilt directions across the vertical projection of the source, strong divergences from the predictions of MBEM (that are likely non-axisymmetric due to the effect of real topography) are found in the region right above the source. We utilize a spatial band-stop filter to remove these numerical effects from the maps in the right column of Fig. 5.

Interesting observations can be deduced from comparison of Figs 4 and 5 on the spatial relationship between (i) source positions leading to the strongest $M_1M/D\_IAMs$ and (ii) areas where the strongest $M_2M/D\_IAM$ occur more often. In particular, deeper sources are more likely to produce important tilt field distortions when their vertical projection to the surface is mostly surrounded by steep topographic features. Indeed, deeper sources that fall below the summit craters and the westernmost edge of the VdB (bottom panel of Fig. 6) induce strong and diffused distortions in correspondence of the slopes to the east and west of the summit zone and of the edge of the VdB (compare the bottom middle maps in Figs 4 and 5; see also the shaded grey area in Fig. 7a). Shallower sources induce more localized effects, mainly driven by the topographic features more closer to the vertical projection of the source itself (compare the top-middle maps in Figs 4 and 5). The above considerations are summarized in Fig. 6. Indeed, moving from bottom (deeper sources) to top (shallower sources) of the colour map (normalized $M_1M/D\_IAMs$ for source positions along the central E–W row of the source grid), the strongest integrated misfits shifts from the centre to the edges of the image, that is from the zone beneath the summit craters to zones beneath highly sloping features of the topography. The upper panel of Fig. 6 also shows how the distortion
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Maps of tilt misfits ($M_2$; see eq. 2) for sets of source positions at different depths. Left column: MBEM/MOGI integrated amplitude misfit ($M_2M/M_IAM$); middle column: MBEM/DVAS integrated amplitude misfit ($M_2M/D_IAM$); right column: numerical/analytical integrated direction misfit ($M_2 IDM$). Numbers reported on each panel represent relative maximum values of the colour scale (the minimum value is always zero; see Table 1 and the text for details).</figcaption></figure>

corresponding to a shallower source is restricted to a smaller area than that corresponding to a deeper source.

From these observations, it follows that topographic perturbations on the tilt field induced by deeper sources will influence more the results of inversion analysis aimed at retrieving the source characteristics using solutions where topography is not properly accounted for. Indeed, the distortion will affect a larger area, implying a more important influence on the input to the inversion scheme.

## 4.3 Influence on the inversion of the source parameters

To explore the above issue, we perform a search of the best source parameters using (i) synthetic data generated through MBEM as the input and (ii) DVAS as the forward model. Relying on the results presented in Figs 4 and 6, we utilize a source position to which important topography-driven tilt field distortions correspond. A source placed below the location of the summit craters, at a depth

of 6 km b.s.l. (test position), fulfills the above condition (see the lowermost panels in Fig. 4 and lower panel in Fig. 6) and also agrees with the findings of previous studies (Aloisi et al. 2011). Fig. 7 shows the tilt field generated by a given pressure increase at this source (whose position is marked by a black star) and calculated through both MBEM (blue arrows) and DVAS (red arrows). As expected, the most important divergences between the two sets of tilt vectors are found in correspondence of the areas of steeper topography around the surface projection of the source (see bottom-middle and bottom-right maps in Fig. 5; see also the shaded grey area in Fig. 7a and curves in Fig. 7b). As the input to the search scheme we utilize data generated by the MBEM at (i) an array of 16 points (array 1) closer to the surface projection of the pressure source and (ii) an array of 16 points (array 2) placed outside the area where the strongest differences between MBEM and DVAS tilt vectors occur (Figs 5 and 7). To exhaustively describe the model space, we search over a 15 × 15 × 13.3 km grid of possible source
<figure><img src="image_8.png" /></figure>

**Figure 6.** Bottom: integrated MBEM/DVAS tilt amplitude misfits ($M_1$; see eq. 1) for source positions along the central E–W row of the source grid (grey dots). Misfit values relative to each source depth are normalized with respect to the average misfit. Misfit values are interpolated via Kriging in the colour map. Note that each point in the colour map represents the normalized rms of the differences between MBEM and DVAS predictions, obtained at each node of the computational grid. Top: MBEM/DVAS tilt amplitude misfits (normalized absolute values) for sources at depths of 1 km (blue curve; source position marked through a blue circle in the bottom panel) and 5 km (red curve; source position marked through a red circle in the bottom panel) below the sea level. The misfits are calculated at the nodes in the central EW row of the computational grid.

positions, with a node spacing of 150 m. The grid is centred on the location of the summit craters and its bottom is placed at 10 km b.s.l.. Nodes above the topographic surface are disregarded. For each *k*th source position in the grid, the source pressure (volume change) is adjusted to minimize the L2-norm of the module of the vector differences between (i) tilt generated through MBEM, with the pressure source in the test position and (ii) tilt calculated through DVAS, with the pressure source at the *k*th grid node. The inverse of the L2-norm values are then converted into the probability density over the whole model space.

Results are shown in Fig. 8. When data generated at the points of array 1 are utilized, the search procedure does not converge to the right solution, due to the important topography-driven distortions that affect most points (Figs 8a and b). Indeed, the most probable solution lies about 3 km from the summit, at a depth of 1 km b.s.l. (white star in Figs 8a and b). Comparison between the adjusted volume change, calculated when the source is in the ‘best’ position, and the volume change assumed when synthetic data were generated through MBEM reveals that the ‘best’ source can account for only 3 per cent of the real volume variation. The real source lies in a low-probability region of the model space (black star Figs 8a and b). When data generated at the points of array 2 are utilized, the search scheme finds a satisfactory solution (white star in Figs 8c and d), whose horizontal position almost coincides with the position of the real source (Fig. 8c) and that is 2 km deeper than the real source (Fig. 8d). This source accounts for 125 per cent of the real volume variation. In this case, the real source lies in a high-probability region of the model space (Figs 8c and d). Note that the maximum source location probability is much lower when data in input to the search scheme are generated at the points of array 1, then when they are generated at the points of array 2.

These findings imply that, under certain conditions, that is pressure source in a ‘critical’ position and observation points in the area where the strongest tilt-field distortions occur, misleading conclu-

<figure><img src="image_9.png" /></figure>

**Figure 7.** (a) Tilt fields due to a spherical source of pressure (black star) at a depth of 6 km b.s.l.. The tilt fields are computed through MBEM (blue arrows) and DVAS (red arrows). The greyscale shading indicates vector difference between tilt fields computed through MBEM and DVAS (light and dark grey represent weak and strong differences, respectively). Points marked in green and yellow belong to arrays 1 and 2, respectively (see text for details). The dashed line encloses points in array 1. (b) E–W cross section passing through the vertical projection of the spherical source of pressure and showing tilt amplitudes (arbitrary units) calculated through MBEM (blue curve), DVAS (solid red curve) and MOGI (dashed red curve).

sions can be drawn on the parameters of the pressure source if tilt data are analysed without accounting for topography.

## 5 CONCLUSIVE REMARKS

We present a strategy to effectively and exhaustively explore the distortions in the tilt field induced by a spherical pressure source, that may arise from topographic features of high-relief volcanoes.

This study is accomplished by comparing the results of a boundary element scheme with analytical solutions and takes, as a case of study, the topographic features of Mt Etna. Nevertheless, the principles that we describe here can be applied to other sites and could involve parameters other than ground tilt, like, for example, horizontal/vertical displacements. We approach the matter in the focus of the present paper in a twofold way. Indeed, we perform extensive screenings of (i) the source positions to which the strongest distortions in the tilt field are associated and (ii) the areas where the strongest distortions in the tilt field are most likely to be observed. Our results show that:

(i) if an analytical solution (MOGI) assuming an elastic half-space with free surface at constant elevation is utilized, the strongest misfit in the amplitude of the tilt vectors are observed around the
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Results of searches for the best source location using tilt data generated through MBEM as the input and DVAS as forward model. In (a) and (b), the input data are generated at the 16 points of array 1 (marked in green; see text for details). In (c) and (d), the input data are generated at the 16 points of array 2 (marked in yellow; see text for details). Black stars: 'real' source. White stars: best inversion source (maximum probability). Black arrows: tilt field induced by the 'real' source of pressure and calculated through MBEM. Red arrows: tilt field induced by the best inversion source of pressure and calculated through DVAS. (b) and (d) present E–W cross-sections passing through the 'real' source. Shading colours indicate source location probability in the whole search domain, with the maximum probability value projected onto the horizontal ((a) and (c)) and vertical ((b) and (d)) planes.</figcaption></figure>

summit of the volcano and when the pressure source is placed below the summit area itself (left columns in Figs 4 and 5);

(ii) if a quasi-analytic method (DVAS) where the vertical distance from the source to the free surface is accounted for is utilized, the strongest misfits in the amplitude of the tilt vectors are observed in correspondence of highly sloping features of the topography and when the source is placed either (1) at greater depths, with its vertical projection surrounded by highly sloping features of the topography, or (2) at shallower depths, with its vertical projection coinciding with highly sloping features of the topography (middle columns in Figs 4 and 5; Fig. 6);

(iii) if either MOGI or DVAS are utilized, the strongest misfits in the direction of the tilt vectors are observed in correspondence of steep slopes and when the source of pressure lies below the uppermost parts of the volcano (right columns in Figs 4 and 5).

We also show that, when the parameter of the deformation source are tentatively retrieved using tilt data, large errors can be introduced if topography is not accounted for. Indeed, if the pressure source lies in a ‘critical’ position (Figs 4 and 6) and the observation points are placed in the area where the strongest tilt misfits are expected (Figs 5 and 7), an inversion procedure based on a half-space analytical formulation, will likely not lead to a satisfactory solution (left panels in Fig. 8). This finding implies that a compromise should

be tentatively established between the use of inversion procedures exploiting numerical 3D-modelling and the need of keeping a reasonable computation time. In this sense, boundary element methods represent an ideal tool (Beauducel & Cornet 1999; Beauducel et al. 2000, 2004; Fukushima et al. 2005).

## ACKNOWLEDGEMENTS

We thank the associate editor, J. Wassermann, and two anonymous reviewers for their constructive comments that helped to improve the initial manuscript. This work was financed through the project MED-SUV that received funding from the European Union’s Seventh Programme for research, technological development and demonstration, under grant agreement No 308665. D. Carbone was partially supported by the IPGP foreigner visiting program. This is IPGP contribution #3602.

## REFERENCES

Aloisi, M., Bonaccorso, A. & Gambino, S., 2006. Imaging composite dike propagation (Etna, 2002 case), *J. geophys. Res.*, **111**(B6), B06404, doi:10.1029/2005JB003908.
Aloisi, M., Mattia, M., Ferlito, C., Palano, M., Bruno, V. & Cannavò, F., 2011. Imaging the multi-level magma reservoir at Mt. Etna volcano (Italy), *Geophys. Res. Lett.*, **38**(16), doi:10.1029/2011GL048488.

Anderson, E., 1936. The dynamics of the formation of cone sheets, ring dykes and caldron subsidences, *Proc. Royal Soc. Edinburgh*, **56**, 128–163.

Anderson, K., Lisowski, M. & Segall, P., 2010. Cyclic ground tilt associated with the 2004–2008 eruption of mount St. Helens, *J. geophys. Res.*, **115**, doi:10.1029/2009JB07102.

Battaglia, M., Gottsmann, J., Carbone, D. & Fernández, J., 2008. 4D volcano gravimetry, *Geophysics*, **73**(6), WA3–WA18.

Beauducel, F. & Cornet, F., 1999. Collection and three-dimensional modeling of GPS and tilt data at Merapi volcano, *Java, J. geophys. Res.*, **104**(B1), 725–736.

Beauducel, F., Cornet, F., Suhanto, E., Duquesnoy, T. & Kasser, M., 2000. Constraints on magma flux from displacements data at Merapi volcano, *Java, Indonesia, J. geophys. Res.*, **105**(B4), doi:10.1029/2009JB07102.

Beauducel, F., De Natale, G., Obrizzo, F. & Pingue, F., 2004. 3-D modeling of Campi Flegrei ground deformations: role of caldera boundary discontinuities, *Pure appl. Geophys.*, **161**(7), 1329–1344.

Calvari, S., Tanner, L.H., Groppelli, G. & Norini, G., 2004. Valle del Bove, eastern flank of Etna volcano: a comprehensive model for the opening of the depression and implications for future hazards, *Geophys. Monogr.*, **143**, 65–75.

Cayol, V. & Cornet, F., 1997. 3-D mixed boundary elements for elastic-static deformation field analysis, *Int. J. Rock Mech. Min. Sci.*, **34**(2), 275–287.

Cayol, V. & Cornet, F., 1998. Effects of topography on the interpretation of the deformation field of prominent volcanoes: application to Etna, *Geophys. Res. Lett.*, **25**(11), 1979–1982.

Cayol, V., Dieterich, J.H., Okamura, A.T. & Miklius, A., 2000. High magma storage rates before the 1983 eruption of Kilauea, Hawaii, *Science*, **288**(5475), 2343–2346.

Charco, M. & del Sastre, P.G., 2014. Efficient inversion of three-dimensional finite element models of volcano deformation, *Geophys. J. Int.*, **196**(3), 1441–1454.

Chaussard, E. & Amelung, F., 2012. Precursory inflation of shallow magma reservoirs at west Sunda volcanoes detected by InSAR, *Geophys. Res. Lett.*, **39**(21), doi:10.1029/2012GL053817.

Cianetti, S., Giunchi, C. & Casarotti, E., 2012. Volcanic deformation and flank instability due to magmatic sources and frictional rheology: the case of Mount Etna, *Geophys. J. Int.*, **191**(3), 939–953.

Crouch, S., 1976. Solution of plane elasticity problems by the displacement discontinuity method. I. Infinite body solution. II. Semi-infinite body solution, *Int. J. Num. Meth. Eng.*, **10**(2), 301–343.

Currenti, G., Del Negro, C. & Ganci, G., 2007. Modelling of ground deformation and gravity fields using finite element method: an application to Etna volcano, *Geophys. J. Int.*, **169**(2), 775–786.

Dieterich, J., Cayol, V. & Okubo, P., 2000. The use of earthquake rate changes as a stress meter at kilauea volcano, *Nature*, **408**(6811), 457–460.

Dvorak, J. & Dzurisin, D., 1997. Volcano geodesy: the search for magma reservoirs and the formation of eruptive vents, *Rev. Geophys.*, **35**, 343–384.

Dzurisin, D., 2003. A comprehensive approach to monitoring volcano deformation as a window on the eruption cycle, *Rev. Geophys.*, **41**(1), doi:10.1029/2001RG000107.

Froger, J.-L., Merle, O. & Briole, P., 2001. Active spreading and regional extension at Mount Etna imaged by SAR interferometry, *Earth planet. Sci. Lett.*, **187**(3), 245–258.

Fukushima, Y., Cayol, V. & Durand, P., 2005. Finding realistic dike models from interferometric synthetic aperture radar data: the February 2000 eruption at Piton de la Fournaise, *J. geophys. Res.*, **110**(B3), doi:10.1029/2004JB003268.

Harrison, J., 1976. Cavity and topographic effects in tilt and strain measurement, *J. geophys. Res.*, **81**(2), 319–328.

Ishii, H. & Takagi, A., 1967. Theoretical study on the Crustal Movements Part. The influence of surface topography (two-dimensional SH-Torque source), The Science Reports of Tohoku University, Fifth Series, Geophysics 19(2), 77–94.

Lisowski, M., 2006. Analytical volcano deformation source models, in *Volcano Deformation*, pp. 279–304, ed. Dzurisin, D., Springer-Verlag.

Lyons, J.J., Waite, G.P., Ichihara, M. & Lees, J.M., 2012. Tilt prior to explosions and the effect of topography on ultra-long-period seismic records at Fuego volcano, Guatemala, *Geophys. Res. Lett.*, **39**(8), L08305, doi:10.1029/2012GL051184.

Masterlark, T., Feigl, K.L., Haney, M., Stone, J., Thurber, C. & Ronchin, E., 2012. Nonlinear estimation of geometric parameters in FEMs of volcano deformation: integrating tomography models and geodetic data for Okmok volcano, Alaska, *J. geophys. Res.*, **117**(B2), doi:10.1029/2011JB008811.

McTigue, D. & Segall, P., 1988. Displacements and tilts from dip-slip faults and magma chambers beneath irregular surface topography, *Geophys. Res. Lett.*, **15**(6), 601–604.

Meertens, C.M. & Wahr, J.M., 1986. Topographic effect on tilt, strain, and displacement measurements, *J. geophys. Res.*, **91**(B14), 14 057–14 062.

Meo, M., Tammaro, U. & Capuano, P., 2008. Influence of topography on ground deformation at Mt. Vesuvius (Italy) by finite element modelling, *Int. J. Nonlinear Mech.*, **43**(3), 178–186.

Mogi, K., 1958. Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them, *Bull. Earthq. Res. Inst. Tokyo*, **36**, 99–134.

Montgomery-Brown, E., Sinnett, D., Poland, M., Segall, P., Orr, T., Zembker, H. & Miklius, A., 2010. Geodetic evidence for en echelon dike emplacement and concurrent slow slip during the june 2007 intrusion and eruption at Kīlauea volcano, Hawaii, *J. geophys. Res.*, **115**(B7), doi:10.1029/2009JB006658.

Ohminato, T., Chouet, B.A., Dawson, P. & Kedar, S., 1998. Waveform inversion of very long period impulsive signals associated with magmatic injection beneath Kilauea Volcano, Hawaii, *J. geophys. Res.*, **103**(B10), 23 839–23 862.

Okada, Y., 1985. Surface deformation due to shear and tensile faults in a half-space, *Bull. seism. Soc. Am.*, **75**(4), 1135–1154.

Peltier, A., Bachèlery, P. & Staudacher, T., 2011. Early detection of large eruptions at Piton de La Fournaise volcano (La Réunion Island): contribution of a distant tiltmeter station, *J. Volc. Geotherm. Res.*, **199**(1), 96–104.

Privitera, E., Bonanno, A., Gresta, S., Nunnari, G. & Puglisi, G., 2012. Triggering mechanisms of static stress on Mount Etna volcano: an application of the boundary element method, *J. Volc. Geotherm. Res.*, **245**, 149–158.

Rizzo, F., 1967. An integral equation approach to boundary value problems of classical elastostatics, *Quart. Appl. Math.*, **25**(1), 83–95.

Segall, P., 2010. *Earthquake and Volcano Deformation*, Princeton Univ. Press.

Trasatti, E., Giunchi, C. & Agostinetti, N.P., 2008. Numerical inversion of deformation caused by pressure sources: application to Mount Etna (Italy), *Geophys. J. Int.*, **172**(2), 873–884.

van Driel, M., Wassermann, J., Nader, M.F., Schuberth, B.S. & Igel, H., 2012. Strain rotation coupling and its implications on the measurement of rotational ground motions, *J. Seismol.*, **16**(4), 657–668.

Voight, B., Hoblitt, R., Clarke, A., Lockhart, A., Miller, A., Lynch, L. & McMahon, J., 1998. Remarkable cyclic ground deformation monitored in real-time on Monserrat, and its use in eruption forecasting, *Geophys. Res. Lett.*, **25**, 3405–3408.

Williams, C. & Wadge, G., 1998. The effects of topography on magma chamber deformation models: application to Mt. Etna and radar interferometry, *Geophys. Res. Lett.*, **25**(10), 1549–1552.

Williams, C.A. & Wadge, G., 2000. An accurate and efficient method for including the effects of topography in three-dimensional elastic models of ground deformation with applications to radar interferometry, *J. geophys. Res.*, **105**(B4), 8103–8120.
SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this paper:

**Figure S1.** Same as Fig. 5, but without applying a spatial bandstop filter to the maps showing the numerical/analytical integrated direction misfit (maps in the right column; see text for details)

(http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggv076
/-/DC1)

Please note: Oxford University Press is not responsible for the
content or functionality of any supporting materials supplied by
the authors. Any queries (other than missing material) should be
directed to the corresponding author for the paper.