Recent Developments and Trends in Volcano Gravimetry

Abstract

The aim of this chapter is to take a look at some developments and new trends in volcano gravimetry. First, we will review the objectives of the research work within this subfield of geophysics, discuss the data and methods it uses, and outline the outputs it strives for. Then, we will turn our attention to three areas where innovative approaches possibly can forward this field of study. The first has to do with the coupling between vertical deformations of the topographic surface (elevation changes) and the observed gravity changes or, in other words, with the removal of the deformation-induced gravimetric signal from the observed gravity changes to obtain the net gravity changes caused by volcanic signals. The second and third areas regard the inversion of the observed gravity changes and deal with two recently or newly developed inversion approaches that both are characterized by the ability to produce a suite of diverse solutions that can be analyzed and discriminated based on additional independent constraints stemming from other earth science disciplines or from the cognition of the interpreter. With this in mind, the final goal is a better understanding of the mechanisms and processes of volcanic unrest or reawakening of a volcano and forecasting the threat of consequent activity and impacts.

1. Introduction

Active and dormant volcanoes cover the face of the earth more plentifully than we usually percept. Their proximity is often richly inhabited by residents or visited by tourists, which gives rise to a tight link between the prosperity, wealth and health, and even the fate of those people and a volcano. As a natural consequence, the knowledge about the possible volcanic threat is of high demand. Our knowledge, however, of volcanic processes and behavior is far from complete or ideal. There is a lot of room for earth science endeavors and research efforts to contribute to this knowledge. Here, we take a look at how much can a specific earth science discipline—gravimetry—contribute to this topic. Gravimetry is a discipline overlapping the fields of geophysics, geodesy, and geodynamics that deals with the geometry, properties, and changes of the earth gravity field by means of its observation, analysis, and interpretation.

When considering volcanic hazards, the knowledge of the following elements is essential: the geological past of a volcano, the geological structure of a volcano and its tectonic setting, and the subsurface processes taking place inside a restless or awakening volcano or deep underneath it within the crust and upper mantle. None of these three elements can be observed directly. The geological past of a volcano is inferred from its products left from the past eruptions and from the morphology of the volcano and its surroundings. The chemical and mechanical composition of the products testifies about the violence, size, and character (such as explosive versus effusive) of past eruptions. Geochemical, petrological, and morphological methods help recover such information. The geological structure of a volcano is not accessible to the naked eye either. Geophysical methods interpreting observations made on the surface, which are sensitive to the distribution of physical properties (parameters) of the volcano edifice, and the crust or uppermost mantle underneath, are applied to compile this knowledge. Seismic methods, gravimetry, magnetics, electrical soundings, induction, self- potential, and ground penetrating radar are examples of such exploration. Similar is the situation with acquiring knowledge on the subsurface volcanic processes linked with migration of magma and associated volatiles. From data monitored or observed on the surface, we aim at reconstructing the reality happening within the volcano. Subsurface physical and chemical (compositional) changes in the magma reservoirs (or changes in the volcano plumbing system in general), in hydrothermal systems, and within their surrounding rock environment, as well as underground mass movement, produce observable geophysical signals on the surface. In the sequel, we shall focus on two such observables: the temporal gravity changes and the vertical displacements (elevation changes) as the vertical component of the deformations of the earth surface.

2. Surface deformations and gravity changes

Processes associated with the physical and chemical changes in the magma storage system and intruding fresh magma both lead to changes in subsurface stress field producing strain that manifests itself on the surface in terms of surface deformations. The subsurface strain field and changes in the temperature field generate changes in the subsurface density distribution, which is sensed on the surface as microgravity changes. Moreover, the movement (migration) of masses such as magma within preexisting feeding or storage systems, as well as formation of new parts of the plumbing system in terms of dikes, sills, inclined or conical sheets, and rings [1], produces another component of changes in the subsurface density distribution. The same holds true for the movement of volcanic fluids (or volatiles in general) within the porous rock medium or along cracks and fractures or for any changes in hydrothermal systems as such. Any changes in the subsurface density distribution manifest themselves on the surface as changes in the gravitational attraction vector. Typically, the vertical component is observed on the surface, referred to as temporal gravity changes. Naturally, both the stress field changes and the subsurface density distribution changes work in a concert, implying a simultaneous occurrence of both the surface displacements and the gravity changes. Therefore, it would be natural and advantageous to interpret surface deformations and gravity changes jointly. We must also mention that these two observables are coupled. This is due to the fact that any surface deformations lead to gravity changes caused by the vertical displacement of the topographic surface in the ambient gravity field regardless of the subsurface mass redistribution.

Volcanological experience has proven that the inflation (or deformation in general) of the volcano edifice is an essential precursor of volcanic eruptions [2–4]. Therefore, the monitoring of surface deformations plays an irreplaceable role in volcanic threat assessment. Surface deformations can be monitored using continuous or repeat survey terrestrial geodetic techniques as well as extraterrestrial techniques [3], such as differential interferometric synthetic aperture radar (DInSAR) or its permanent scatterer modification (PSInSAR). Gravity changes turned out to be also a valuable tool for predicting the reactivation of a volcano or for studying volcanic unrest [5, 6]. Cases have been reported when gravity changes were detected before seismic or other precursors at a reawakening volcano (e.g., [2, 7–9]). Gravity changes can be monitored continuously using tidal relative gravimeters or absolute gravimeters or in a repeat campaign (survey) mode at a network of gravity points (benchmarks, stations) covering a representative area over the volcano. While the repeat campaign (time-lapse) observations can better capture the spatial extent of the gravity changes, the continuous measurements at a station (or several stations) can pick up the fast changes to which the repeat survey mode is blind in the intracampaign period. The best, of course, appears to be the combination of the two modes [5]. Here, we shall focus on the time-lapse gravity changes and their interpretation. The continuous gravity changes and the benefits of their interpretation are well described in [5, 6].

Surface deformations are sensitive particularly to source pressure or volume changes [3, 4, 10]. Gravity changes help characterizing the source process by sensing mass changes along with volume changes, thus pointing at the source density changes [11–17]. Monitoring and interpreting temporal changes in 3D gravity data has been recently referred to as “4D (micro)gravity” or 4D gravimetry [5, 6]. A comprehensive overview of the purpose and benefits of the 4D gravity, both the discrete (time-lapse, repeat survey) and the continuous, illustrated on several cases at specific volcanoes, is given in [5]. A special case happens when gravity changes are observed without any significant surface deformations or vice versa [12, 15, 18]. By “significant” we mean within the accuracy threshold of the measuring technique. This strengthens the argument for concurrent monitoring of both the surface deformations and the gravity changes in volcanic areas.

Surface deformations and gravity changes can be interpreted either as stand-alone quantities or jointly. We have already indicated the pitfalls of monitoring and interpreting just one of them and the advantages of monitoring and interpreting them jointly. The (sole or joint) interpretation is based on either modeling or inversion. Modeling is based on assuming a hypothesis about a source and its parameters, simplifying the structural rock environment, and solving the forward problem yielding the model surface deformations and model gravity changes. By iteratively changing the source parameters or the hypothesized model source, modeling is carried out until a reasonable fit (“match”) is achieved between observed and modeled quantities (displacements and gravity changes). When a fit is reached, the model source with its parameters represents one of the possible solutions representing the volcanic process taking place at depth. Two approaches to modeling are at hand: analytical and numerical.

In analytical modeling, significant simplifications are applied to the source as such (its geometry) and to the structural environment: neglecting the topography and subsurface geologic structure (assuming the environment to be a homogenous isotropic elastic or viscoelastic half-space). Presuming simple geometries of the sources (such as magma chambers, dikes, and sills) leads to closed-form analytical expressions linking source parameters and the observable surficial quantities. Hence, analytical modeling can give quick results and first impressions about the nature and characteristics of the source process. On the contrary, the oversimplifications used in analytical modeling can lead to severe distortions of the obtained results with the risk of misleading the interpreter. A handy overview of analytical modeling approaches for pressurized source geometries given as spherical, ellipsoidal, cylindrical, sill-like, and dike-like are given in [6].

In numerical modeling, the simplifications of the analytical approach are remedied. The numerical approach takes into account the effect of topography, structural discontinuities and inhomogeneity, and even rheology. Such approach calls for applying finite element or boundary element methods (FEM, BEM). An instructive overview of the numerical modeling techniques in volcano geodesy is given again in [6].

Inversion works differently: a nonlinear inverse problem is solved by means of computing the source parameters directly from the observed surficial data. This gravimetric inverse problem is nonunique and ill-posed. The methods to solve it are the same as those in potential field inverse problems.

Some volcanic areas, especially calderas, display during their restlessness a special behavior—they produce systematic trends in their gravity/height ratios. By analogy, these ratios are typically referred to as gravity/height gradients. When the plots of gravity changes versus elevation changes show systematic trends clearly distinguished by two linear boundaries [free air gradients (FAGs) and Bouguer gradients], they are considered as signatures of the dominating typified subsurface processes [7, 8, 19–24]: magma rejuvenation or drainage, vesiculation or degassing, and water table rise or fall.

The mathematical source solutions (model solutions) obtained from either modeling or direct inversion must be associated with processes taking place in the volcano edifice or even deeper in the crust or lithospheric mantle and furthermore often with a volcanic threat or a prediction about an eruption. This is not a straightforward task. However, this topic is outside of the scope of our chapter here. We shall pay a particular attention only to inverting and interpreting temporal gravity changes.

There is an everlasting problem present with the monitoring and interpretation of geodetic (deformations) and gravimetric (gravity changes) data in volcanic areas: they are neither easy nor cheap to observe—their temporal as well as spatial resolution is never high enough to satisfy the inversion and interpretation demands. This will become particularly evident when we shall further below discuss the novel approaches to gravity inversion that require the input data given on an equidistant (regular) and dense enough grid. When the spatiotemporal gravity changes do not have a sufficient resolution and accuracy, it becomes hard, nay impossible, to discriminate among the possible sources of volcanic unrest and possibly draw reasonable conclusions about the threat of impending eruption.

3. Decomposition of superimposed gravity signals

Several natural (physical) phenomena produce a change in gravity observed at a point (benchmark) on the earth surface: changes in atmospheric attraction (due to mainly atmospheric pressure changes), tidal effects such as solid earth tides and ocean (loading) tides, hydrological effects such as groundwater table level changes, changes in snow or ice cover, soil moisture, volcanic effects associated with changes in magma storage (state of magma and magma transport), as well as magmatic volatiles/fluids. All of these signals (components)—gravitational effects of the individual contributions—are superimposed to form the gross observed gravity changes. When we aim at studying a particular geodynamic phenomenon, such as magma replenishment into an existing plumbing system of a volcano, a magma intrusion, or eventually magma vesiculation/degassing in a shallow chamber or vent, using gravimetry, we must be able to decompose the composite gravity signal. In other words, we need to remove (correct for) all unwanted signals. The net signal, stripped of any unwanted components, is typically referred to as “residual gravity changes”.

First, we remove the contribution of the environment: gravity change components imposed by the atmosphere and tides. Then, we need to deal with the hydrological component due to pluvial water (precipitation). This turns out to be a cumbersome correction to handle in practice, the size of which often is at the level of the net signal chased after. Magma is associated with its volatiles and brines. Hydrothermal systems often reside above the magma storage systems. Physical changes (such as temperature and pressure) in both the magmatic and the hydrothermal systems have gravimetric signatures that may be difficult, nay impossible, to discern. Similar situation may arise with the gravimetric signatures associated with the transport (mobility) of magma and hydrothermal fluids.

This has become obvious in monitoring and interpreting unrest at calderas. Sometimes, researchers come to diverse conclusions interpreting the same data (gravity changes and surface deformations). Some interpretations favor a magma intrusion process (e.g., [2, 19, 25, 26]), whereas the opposing ones prefer a hydrothermal reasoning (e.g., [11, 14, 27–29]). Naturally, these two source processes may act simultaneously complementing each other, being referred to as hybrid unrest (e.g., [15, 30, 31]).

For the purpose of inverting and interpreting the net signal, the residual gravity changes can be compiled using the following decomposition [32] (see also Figure 1):

$\text{Δ}{g}^{res}=\text{Δ}{g}^{obs}-\text{Δ}{g}^{ext}-\text{Δ}{g}^{inst}-\text{Δ}{g}^w-\text{Δ}{g}^{def}-\text{Δ}{g}^{surf}$ E1

where the superscripts denote the following components of the gravity change: residual (res), observed (obs), external (ext), instrumental (inst), hydrological (w), deformation induced (def), and surficial (surf).

By external component, we mean environmental effects composed of tidal and atmospheric effects. Tidal effects comprise solid earth (body) tides and ocean loading effects. Atmospheric effects are gravity changes induced by pressure and temperature changes in the atmosphere (atmospheric attraction and loading effects). These effects are known and respective corrections have been published (e.g., [33–37]).

As instrument/survey effects, we consider the drift of gravimeters and adjustment of redundant measurements. We want to highlight one particular effect that must not be neglected in monitoring and interpreting time-lapse gravity change observations. Various gravimeters have various heights of the sensor, measured from the bottom of the instrument. In addition, various plates or tripods can be used, implying various heights of the bottom of the instrument above a benchmark (gravity station, gravity point). If various types of gravimeters or plates/tripods are used within an epoch survey or between epochs, observed gravity must be reduced to the elevation of the benchmark. Otherwise, significant systematic errors could be introduced to the residual gravity changes. To reduce the meter reading to the ground, actual vertical gradient of gravity (VGG) is needed. It may be either observed at each station, or estimated, as will be discussed further below (cf. also [38]).

Under hydrological effects, we mean the changes in the groundwater table level as well as in the soil/rock moisture (due to precipitation or drought). Strictly speaking to correct for this environmental signal, detailed 3D modeling needs to be performed, which is typically unachievable due to lack of required input data and/or knowledge about the near surface (or subsurface) geological structure (porosity, impermeable interfaces, etc.). Various approximations and estimates are used in practice; data on water table from wells are used wherever available. The simplest approximation is a planar Bouguer effect of a water table level change (e.g., Appendix in [14]). The hydrological correction is not the subject of our interest here, so we refer the reader to published works (e.g., [39–41]).

Deformation-induced topographic effects (DITEs) are of our high concern here. They are treated in Section 5.

Figure 1.

Schematic sketch of the sources contributing to the composite gravity changes.

Under the surface effects, we mean gravitational effects of mass changes taking place on the surface, but unrelated to surface deformation, such as magma extrusion, dome growth, dome collapse, lava flow, lahar, and flank collapse. These masses can be surveyed on the surface, their shapes and sizes can be digitized, and their gravitational effects can be computed by a numerical realization of the Newton volume integral for the vertical component of the attraction and subsequently subtracted (corrected for). This may be achieved with the help of photogrammetry, laser scanning, LIDAR, etc. These corrections are out of scope of our interest here. We refer the reader to published works (e.g., [16, 42]).

The residual gravity changes are consequently related to the subsurface volcanic processes associated with magma migration (rejuvenation or drainage of existing reservoirs and feeders, intrusions, and propagation along new paths, such as diking), hydrothermal fluid (volatile) migration, density changes due to physical and chemical (compositional) changes of the residing magma (cooling, heating, vesiculation, degassing, fractional crystallization, differentiation, mixing, partial melting, etc.), and density changes of the geological structure of the edifice due to stress-induced strain field (including the deformation of density interfaces) and due to temperature field changes via thermal expansivity. All of these changes represent the net sought signal, which in its nature is a fairly complex composite signal itself. This component of the gravity changes, the “residual gravity changes”, is the subject of inversion or forward modeling and of subsequent interpretation of magma-related processes.

For decompositions analogical to our Eq. 1, see, e.g., Eq. 1 and Figure 1 in [6], Eq. 1 and Figure 1 in [16], and also [13, 19, 21, 43–47].

4. Recent developments and new trends

We want to take a look at some of the recent developments and contributions to the field of volcano gravimetry in the sequel. We will turn our attention to three topics:

1. Revisiting the DITEs on gravity changes dealing with the coupling between vertical displacements (elevation changes) and gravity changes, attempting to propose a more rigorous and accurate way of handling them (Section 5);

2. Taking a look at the benefits arising from applying a novel inversion methodology, as applied to time-lapse gravity changes, to which we shall briefly refer as the “Prutkin inversion methodology” (Section 6); and

3. Probing the possibilities of applying a novice inversion methodology, developed by Pohánka, based on n-harmonic (polyharmonic) functions, referred to as the “harmonic inversion method”, to time-lapse gravity changes (Section 7).

For all the three above topics, we chose the Central Volcanic Complex (CVC) of Tenerife, Canary Islands, as our case study playground (Figure 2), on which we ran either synthetic simulations based on high-resolution digital elevation models and Mogi point sources or processed and inverted the gravity changes observed during the 2004 to 2005 volcanic unrest. The selected area has a suitable (significant and jaggy) topography for our simulations, as it comprises a caldera at the average elevation of roughly 2000 m asl, within which twin stratovolcanoes Teide and Pico Viejo are located reaching altitudes 3718 and 3135 m asl, respectively. Among other unrest indicators, spatiotemporal gravity changes were observed at the CVC at 14 benchmarks of a rapid reaction network between May 2004 and July 2005. No statistically significant areal surface deformation (either inflation or deflation) was observed accompanying these gravity changes. The observed gravity changes were corrected for tidal and hydrological effects. These point gravity data were taken as input data in both our inversion approaches (Sections 6 and 7). They were interpolated (extrapolated) onto a regular (equidistant) grid.

Figure 2.

(a) Location of the Canary Islands and Tenerife and (b) the caldera and the twin stratovolcanoes Teide and Pico Viejo (PV).

5. Coupling between vertical deformations and gravity changes

Let us consider the gravitational effects imposed by the deformation of the topographic surface, i.e., by vertical displacements (elevation changes), illustrated (Figure 3), for instance, on inflation (uplift).

Figure 3.

Schematic sketch of the DITEs on the gravity change. Benchmarks 5 to 15, unlike 1 to 4, experienced a significant uplift (Δ_hi). “A” denotes the predeformation topographic surface, whereas “B” denotes the postdeformation surface.

During the deformation (uplift), the gravity station (benchmark) is vertically displaced together with the topographic surface. We note that the subsurface density (geological) structure (including discontinuity interfaces) is deformed, too, but we treat this effect as part of the residual gravity changes to be inverted and interpreted. The reason for this is that, unlike the topographic surface, which is known (described by a DEM), the geological structure is usually unknown or poorly known. We can hypothetically split the ground deformation effect into two subsequent stages (steps) with their respective effect components. In the first step, we move (lift up) the benchmark vertically from its original position (P) on the predeformation topographic surface to its new postdeformation position (P*) in the ambient gravity field (in “free air”) without deforming (inflating) the topographic surface (topographic masses) yet. In the second step, we move (deform, inflate) the topographic masses. Thus, we can write for these two separated components:

$\text{Δ}{g}^{\text{def}}\left(P^{\text{*}}\right)=\text{Δ}{g}^{\text{FAE}}\left(P^{\text{*}}\right)+\text{Δ}{g}^{\text{TDE}}\left(P^{\text{*}}\right)$ E2

where $\text{Δ}{g}^{\text{def}}$ is the DITE; $\text{Δ}{g}^{\text{FAE}}$ is the gravity change due to the vertical displacement of the benchmark in ambient gravity field (in “free air”), in the sequel referred to as the “free air effect” (FAE); and $\text{Δ}{g}^{\text{TDE}}$ is the gravitational effect (vertical component of the attraction vector) of the topographic masses that were displaced into the volumetric domain between the predeformation and the postdeformation topographic surfaces (green domain in our schematic sketch of Figure 3). We have coined this effect the “topographic deformation effect” (TDE). Already from the sketch of Figure 3, we can intuitively deduce that, at benchmarks 1 to 4, which experienced no uplift (no significant elevation changes), the FAE will be zero, whereas the TDE will not be zero. The masses displaced into the topographic deformation rind (the green domain in our sketch) will impose gravitational effect not only on benchmarks 5 to 15, which have undergone nonzero elevation changes, but also on benchmarks 1 to 4, which experienced no uplift, zero elevation changes. This is an important insight. The size of the TDE at a specific benchmark will depend on its vertical position relative to, and horizontal distance from, the deformation shell (rind) and on the volume and shape of the deformation shell. Below, we shall examine the TDE by numerical simulations. As a consequence, the DITE might be nonzero even at benchmarks with zero elevation changes. This may hint on the possible insufficiency of expressions of the DITE that linearly depend on the vertical displacement, i.e., those being a product of some sort of a gradient (such as Bouguer or free air) and the vertical displacement. These issues associated with the DITE will be examined by numerical simulations further below.

5.1. FAE and its approximations

The FAE amounts to the vertical displacement of the benchmark ( $\text{Δ}h\left(P\right)$ ) times the true (in situ) VGG:

$\text{Δ}{g}^{\text{FAE}}\left(P^{\text{*}}\right)={\left(\partial g/\partial h\right)}^o\left(P\right)\text{*Δ}h\left(P\right)$ E3

Inevitably, the VGG must be observed in situ at the benchmark (hence the superscript “o”). This can be practically achieved by relative gravimeters observing in a so-called tower mode, i.e., on the ground and at a certain height above the benchmark, such as 1 m, on a tripod (e.g., [38]).

If the VGG is not observed in situ, and its value remains unknown, it is usually (one could say habitually) approximated by the “theoretical FAG”, also called the “normal FAG” (e.g., [43, 48]), which is the constant average vertical gradient of normal gravity at the surface of the normal reference ellipsoid. In the sequel, we shall abbreviate the normal FAG as NFAG (NFAG=-308.6 μGal/m). Observations indicate that the approximation of the true VGG by the NFAG can introduce a relative error of up to 88% in rugged terrain of alpine mountains, such as the High Tatras of Slovakia [38], or 77% at the CVC of Tenerife as shown below.

Figure 4.

TNFAG values at the CVC of Tenerife. Minima (in absolute value) correspond with sharp concave topofeatures (valleys, gullies, drains), whereas maxima (in absolute sense) with sharp convex topofeatures (peaks, ridges, pillars, crests, and ribs). Max = −70, min = −481, mean = −315, SD = 34 (μGal/m).

Therefore, we propose a better approximation of the true VGG, in case it is not observed, than the NFAG [32]. In rugged terrain regions, the strongest (primary) contribution to the VGG comes from the terrain (topography and eventually also bathymetry for near shore points) from within the near vicinity of the benchmark. This contribution can be modeled (adopting a constant topodensity) using precise high-resolution digital terrain models (ibid). The NFAG can be refined (corrected for) using this topographic contribution to the VGG. In the sequel, we shall abbreviate the “topocorrected NFAG” as “TNFAG”. We show numerical values of the TNFAG for the CVC of Tenerife in Figure 4. The computational method and the used DEM are described in [32]. The topographic contribution (correction) amounts in the CVC area of Tenerife to 77% (in relative sense) of the NFAG itself, which is highly significant. This is the amount by which the true VGG may differ from the NFAG in a volcanic area of rugged topography.

The secondary contribution to the true VGG comes from the underground geological structure (density anomalies) of the earth. This contribution typically remains unmodeled, whereas the structure remains unknown. Alternatively, Rymer [21] and Berrino et al. [43] propose an approximation (of the VGG), a refinement to NFAG, based on the Bouguer anomaly map of the area. The contribution of geology to FAG in volcanic areas with significant (rugged) topography is expected to be smaller than that of topography, although, in flat areas (such as some calderas), the situation may be vice versa.

5.2. TDE and its approximations

The TDE is the gravitational effect (attraction) of the masses between the topographic surfaces before and after the deformation (within the topographic deformation shell) on the observed gravity change at each benchmark of the survey (see Figure 3). It must be computed by numerical evaluation of the Newton integral (for vertical component of the attraction vector) over this volumetric domain with assumed mean constant topographic density ( ${\rho }_0$ representing the near surface rock environment).

The TDE is in volcanogravimetric applications commonly approximated by a Bouguer plate effect (e.g., [21, 43]), referred to as planar “Bouguer deformation effect” (BDE): ${\text{Δg}}^{\text{TDE}}({\text{P}}^{\text{*}})\approx {\text{Δg}}^{\text{BDE}}({\text{P}}^{\text{*}})$ , where ${\text{Δg}}^{\text{BDE}}\left({\text{P}}^{\text{*}}\right)=2\text{πG}{\rho }_0\text{*Δh}\left(\text{P}\right)$ G being the gravitational constant. Alternatively, the TDE is approximated by a spherical Bouguer term (e.g., [21, 49]), which is a combination of three dilation effects respective to a Mogi point source model [32]: ${\text{Δg}}^{\text{TDE}}({\text{P}}^{\text{*}})\approx (4/3)\text{πG}{\rho }_0\text{*Δh}\left(\text{P}\right)$ . For a detailed discussion about these approximations, as well as how the Bouguer gradients respective to these two Bouguer approximations relate to the gravity/height gradients, we refer the reader to [32]. In both (planar and spherical) Bouguer approximations of the TDE, the BDE depends linearly on the displacement of the benchmark. This fact alone is a clear indication that the BDE might not be a sufficient approximation of the TDE: for benchmarks with zero vertical displacement (such as benchmarks 1–4 in Figure 3), yet close enough to the deformed topographic masses to sense their attraction (at the order of 10 μGal), the BDE will be zero, whereas the TDE will not. This difference between the TDE and its BDE approximation (combined with the departure of the in situ VGG from the NFAG) explains the occurrence of observed data at many volcanoes that in terms of gravity/height plots often dramatically and erratically deviate from the linear trends of the FAG or Bouguer-corrected FAG (BCFAG).

To illustrate the size and spatial characteristics of the TDE—most of all that it significantly departs from its Bouguer approximation (BDE), we computed the TDE at the CVC of Tenerife for a synthetic deformation (surface displacement) field (shown in Figure 4 of [32]) generated by two (shallow and deep) Mogi sources and compared it to the planar BDE for the same displacement field by displaying the difference between the two in Figure 5. The shallow Mogi source was located at the depth of 500 m roughly below the Teide summit, scaled to have a displacement magnitude of 1 m, whereas the deep source was located at the depth of 6 km roughly 5 km to the northwest of the twin stratovolcanoes, scaled to have a displacement magnitude of 50 cm. For such significant surface deformations, the difference between the TDE and the BDE is also significant (see Figure 5 and statistics given in the captions to the figure). In future simulations, we plan to compute the TDE also for smaller deformation fields to assess its deviation from the BDE approximation for smaller displacements.

Figure 5.

Difference between the TDE and its approximation by the planar BDE, both respective to a synthetic displacement field simulated at the CVC of Tenerife. Max = 56, min = −60, mean = −1, SD = 6 (μGal).

5.3. DITE and its approximations

The DITE on gravity change is the sum of the FAE and the TDE. Because we did not have observed values of the VGG available for our study area, we used the TNFAG values as an “allegedly sufficient” approximation of the VGG in the first (FAE) term of DITE. The second (TDE) term was computed numerically via the Newtonian integration. Next, we compared the DITE computed in this fashion (TN-DITE strictly speaking) to two of its commonly used approximations [32]: (a) the NFAG approximation of DITE ( ${\text{Δg}}^{\text{def}}\left({\text{P}}^{\text{*}}\right)\approx \text{NFAG*Δh}\left(\text{P}\right)$ ) and (b) the BCFAG approximation of DITE ( ${\text{Δg}}^{\text{def}}\left({\text{P}}^{\text{*}}\right)\approx \text{BCFAG}\left({\rho }_0\right)\text{*Δh}\left(\text{P}\right)$ ). This comparison leads to interesting findings: (1) the sharp short-wavelength spatial variability in the FAE term caused by the contribution of topographic masses to the VGG (the difference between NFAG and TNFAG) cancels out to some (if not great) extent with its counterpart in the TDE term (the difference between TDE and BDE) and (2) the NFAG approximation of DITE works better for narrow sharp displacement fields (generated by shallow point-like sources), whereas the BCFAG approximation of DITE works better for displacement fields of larger horizontal extents of several kilometers (generated by deeper sources). In other words, the situation with local (sharp/narrow) deformation is closer to ”normal free air” situation, whereas the larger (more widespread) areal deformation is better approximated using the Bouguer plate effect in addition to the normal free air (“Bouguer situation”).

6. Prutkin inversion methodology

The methodology under a brief working title “Prutkin inversion methodology” refers to a modular inversion methodology that consists of several modules or, in other words, several processing and inversion steps or algorithms. This modularity leads to producing several diverse inversion solutions that all equally well satisfy (mathematically) the input data (time-lapse gravity changes in the case of volcano gravimetry). The methodology was originally developed and applied to structural geophysical studies at global, regional, and even local scales [50–53]. We [31, 54] have applied it recently to temporal gravity changes and the results seem promising.

This methodology is used to be referred to as the “method of local corrections”, which is not telling, as the name refers to the last algorithm of the methodology only. That is why we refer to it here loosely as the “Prutkin methodology”. The methodology consists of the following blocks or algorithms: (a) trend removal (optional), where trend is defined as a 2D harmonic function and computed from the input data; (b) decomposition into shallow/deep (or more depth interval) field components (optional) using a triple harmonic continuation numerical procedure; (c) 3D line segments approximation of the sources, where the gravitational effect of these sources approximates the input gravity data; (d) field decomposition based on the individual or grouped line segment effects; and (e) solutions in terms of potato-shaped (closed compact star-convex) source bodies and/or contact surfaces (density discontinuity interfaces). All of these data processing and inversion steps are described mathematically in [31, 53]. When the methodology is applied to (time-lapse) gravity changes instead of gravity anomalies, the represent subsurface time-lapse (temporal) density changes instead of density anomalies. For instance, the potato-shaped bodies are volumetric domains of a homogenous density change.

A pilot study was carried out on gravity changes observed at Mayon Volcano, Philippines [54]. This study revealed that the methodology is applicable to time-lapse gravity changes observed during unrest or volcanic activity. The second study [31] was performed on gravity changes observed during the 2004/5 unrest on Tenerife [15]. This study manifested the true potential of the methodology. It produced a more revealing gravimetric picture of the unrest compared to the previous gravimetric interpretations. The decomposition of the observed gravity changes into shallow and deep fields lead to the unveiling of the deeper source that was interpreted as magmatic [31]. This was possible only thanks to the unique feature of the methodology—the depth-wise decomposition of the signal based on the triple harmonic continuation. Previous interpretations were not capable of detecting this deeper source characterized as of magmatic origin. The sources pertinent to the shallow field were interpreted as hydrothermal fluids [31]. Consequently, the 2004/5 unrest on Tenerife was evaluated, based on the gravimetric interpretation, as hybrid (due to the presence of both magmatic and hydrothermal sources). Seismic constraints were adopted to assist and strengthen the gravimetric interpretation [31].

The great advantages of the Prutkin inversion methodology dwell in its modularity and versatility. Several combinations can be put together in processing the input data. Several combinations in decomposing the input signal can be introduced. This leans on whether trend removal is performed or not, whether the input signal is decomposed into shallow and deep fields (or even several depth components) or not, and what depth level is chosen for its decomposition. If the input data allow the approximation of sources by several line segments, additional decomposition of the input signal (field) is possible by separating the signal respective to the individual line segments or various groups of line segments. In addition, each field component can be inverted either in terms of a density contrast contact surface (interface) or in terms of 3D star convex compact homogenous (“potato-shaped”) source body. Moreover, a priori density contrasts (differential densities), as well as a priori depths for the asymptotic planes of the interfaces, have to be selected, as they are free parameters. The choice of various density contrasts of various source bodies and/or interfaces creates additional variability in the sizes and shapes of the inversion solution sources. Consequently, the versatility and modularity of the methodology results in a possibly rich set of diverse inversion solutions. All of these solutions will equally well satisfy the input data (in terms of the forward problem). However, not all of them will be equally meaningful or reasonable in terms of geological, tectonic, geodynamic, or volcanological considerations. The great advantage of the presented methodology dwells in its ability to produce plentiful diverse solutions to the gravimetric inverse problem that can be at a later stage discriminated based on independent geophysical constraints (such as seismic) or cognition stemming from other earth science disciplines (such as geochemistry and petrology) or alternatively stemming from the experience of the interpreter. An illustration of such variability of gravimetrically acceptable diverse solutions can be found for the Kolarovo gravity anomaly in the Danube Basin presented in [53]. The same applies to gravimetric inversion of time-lapse gravity changes.

7. Harmonic inversion method

The harmonic inversion method was developed by Pohánka [55, 56]. It has been applied to gravity data in structural geophysical studies. We were curious about its applicability to temporal (time-lapse) gravity changes, so we stepped to testing it on the gravity changes observed during the 2004/5 unrest on Tenerife [57]. We used the same data set as in the Prutkin inversion methodology (Section 6) to be able to compare inversion results of both methodologies.

The gravimetric inverse problem is nonunique. It has infinitely many solutions. Some special solutions can be sought, characterized by being smooth. This methodology searches for subsurface 3D solutions satisfying three conditions [55, 56]: (1) the solution is a n-harmonic function, where n is a small integer (n=3, n=4); (2) the solution has an extremum preservation property for a point source field, i.e., the main extremum of the solution coincides with the location of the point source generating the input gravitational field; and (3) the solution can be expressed as a linear integral transformation of the input gravity field. The tetraharmonic (n=4) solutions are called “characteristic density” or shortly “chi density”. Such harmonic solutions are smooth and unrealistic, as they do not reflect subsurface geological settings. Therefore, in the subsequent step, these harmonic solutions are transformed to piece-wise homogenous density distributions by an automated iterative procedure, in which the harmonic solution (as a 3D subsurface function computed from the 2D surficial input data) serves as a mediator controlling the iterations in terms of the goodness of fit. The piece-wise homogenous solution consists of a set of homogenous closed-form source bodies of specific density contrasts with respect to an ambient constant density. Each source body can have its own density contrast.

The current version of the methodology works with a triharmonic function (n=3) called “quasi-gravitation”, which serves the needs of the mediator in the iterative process resulting in the piece-wise constant subsurface density distribution, given as a set of homogenous source bodies. Density contrasts (differential densities) of these source bodies have to be selected a priori by the interpreter before running the iterative procedure. The iterative procedure, which modifies the shapes and sizes of the individual source bodies, runs until a satisfactory (preselected sufficient) fit is obtained between the quasi-gravitation of the found source bodies and the quasi-gravitation of the observed input gravity data. This fit guarantees also the fit between the surface gravity of the found source bodies and the observed input gravity. Regarding all the details about this inversion methodology based on n-harmonic functions, we refer the reader to [57]. When inverting gravity, the differential densities of the solution represent structural density anomalies. When inverting temporal (time-lapse) gravity changes, the differential densities represent time-lapse density changes.

The residual gravity changes of the 2004/5 Tenerife unrest (14 point data) were extrapolated onto a regular (equidistant) grid 60×60 km with a step of 200 m [57]. Quasi-gravitation (as a subsurface 3D function) was computed from the gravity changes in a domain 40×40 km to the depth of 12 km (2,424,060 cubic cells of the 200 m size). The computed quasi-gravitation had 11 local extrema (7 positive and 4 negative). These extrema served as seeds for growing the homogenous source bodies during the automated iterative process. For each extremum/body, the interpreter has to preselect the value of the differential density. The sizes and shapes of the obtained source bodies depend not only on the observed input data but partly also on the selection of the individual differential densities. Therefore, several different solutions, in terms of a set of source bodies (equal in number to the number of the extrema of the quasi-gravitation), can be obtained for various selections of the differential density of each body. In our pilot study [57], in which we tested the applicability of the methodology to gravity changes, we focused on investigating the variability of the solutions depending on the choices of the individual differential densities. Several solutions originating from the harmonic inversion approach were presented for the gravity changes of the 2004/5 Tenerife unrest in [57]. We have learned the following lessons from the pilot study.

Figure 6.

Harmonic inversion solution for the 2004/5 Teide unrest: source bodies with differential density 1 kg/m³.

The smaller the differential density of a specific source body is, the larger is the source body and the more complex is its shape. The larger the differential density of a body is, the smaller and more spherical is the source body (tending towards a point mass). Interesting shapes and connections among the source bodies possibly indicating pathways of higher permeability were achieved for differential densities as low as 1 kg/m³ (cf. Figure 6).

One of the admissible solutions obtained by the inversion approach based on harmonic functions is illustrated in Figure 6. For such small differential densities, the sizes of the source bodies might be unrealistically overgrown. We have yet to learn how to interpret (volcanologically) the source bodies found by the harmonic inversion approach. Some hints can be already grasped from our pilot study. One option is to view them as volumetric zones delineating the possible presence of magmatic sources of higher differential densities and lesser dimensions discretely distributed within the zones (such as dike/sill/sheet swarms or magma batches [1]) and/or zones of density changes due to migration of hydrothermal fluids through porous media. A follow-up work is needed to pursue this issue.

The case study based on the 2004/5 time-lapse gravity changes of the CVC on Tenerife indicates that the harmonic inversion methodology appears promising for inverting and interpreting gravity changes observed in volcanic areas during unrest or reactivation. In the follow-up work, we shall attempt to assign volcanological interpretation to the source bodies obtained by this inversion approach, such as the one presented in Figure 6. Additional case studies, in general, are needed to establish the link between the solutions obtained by the harmonic inversion method and their physical interpretation in terms of processes taking place inside the volcano edifice or deeper below.

8. Discussion and conclusions

The Achilles heel of the two presented inversion methodologies is that they need input data at an equidistant grid and, moreover, a fairly fine step to produce sources with reasonably detailed geometry. The data available in volcano gravimetry are typically sparse point data. Consequently, interpolation/extrapolation is inevitable. However, the extrapolation can easily introduce false signal features into the input, or artificially exaggerate the amplitude of highs and lows in the data, distorting the picture about the sources. There is no way around other than densifying the networks and acquiring more observations during the campaigns.

Both the inversion methods dealt with here require that the input data be given on the surface of a half-space (on a plane), which was not the case of the Tenerife gravity changes observed on and referred to the topographic surface. In our case studies reported here, this issue was neglected and the data were treated as if given on a plane at a level of roughly 2000 m asl. Strictly speaking, the data should have been harmonically continued from topography to a level surface, which in practice is numerically hardly achievable for sparse point data.

In potential field inverse problems, we often deal with composite signals, whereas the decomposition into multiple sources inherently remains an ambiguous task. Due to the nonuniqueness of the inverse gravimetric problem, also the solutions found by both the Prutkin inversion methodology and the harmonic inversion are again among many admissible solutions. Additional constraints from other earth science disciplines are required to validate physical (geological, volcanological) feasibility of any such solution. However, it is of great advantage when gravimetry alone provides several sets of admissible solutions that later on may be discriminated using available geologic or geoscientific constraints in the search for the most probable realistic scenario.

Acknowledgments

This work was supported by the Slovak Research and Development Agency (APVV) under Project No. APVV-0724-11 as well as by the VEGA grant agency under Project Nos. 2/0042/15 and 1/0141/15.