Open AccessReview

Modeling 3D Free-geometry Volumetric Sources Associated to Geological and Anthropogenic Hazards from Space and Terrestrial Geodetic Data

Antonio G. Camacho

and

José Fernández

Institute of Geosciences (CSIC-UCM), C/ Doctor Severo Ochoa, 7, Facultad de Medicina (Edificio entrepabellones 7 y 8, 4a planta), Ciudad Universitaria, 28040 Madrid, Spain

Author to whom correspondence should be addressed.

Remote Sens. 2019, 11(17), 2042; https://doi.org/10.3390/rs11172042

Submission received: 26 July 2019 / Revised: 22 August 2019 / Accepted: 22 August 2019 / Published: 29 August 2019

(This article belongs to the Special Issue New Frontiers of Multiscale Monitoring, Analysis and Modeling of Environmental Systems via Integration of Remote and Proximal Sensing Measurements)

Download

Browse Figures

Versions Notes

Abstract

Recent decades have shown an explosion in the quantity and quality of geodetic data, mainly space-based geodetic data, that are being applied to geological and anthropogenic hazards. This has produced the need for new approaches for analyzing, modeling and interpreting these geodetic data. Typically, modeling of deformation and gravity changes follows an inverse approach using analytical or numerical solutions, where normally regular geometries (point sources, disks, prolate or oblate spheroids, etc.) are assumed at the initial stages and the inversion is carried out in a linear context. Here we review an original methodology for the simultaneous, nonlinear inversion of gravity changes and/or surface deformation (measured with different techniques) to determine 3D (three-dimensional) bodies, without any a priori assumption about their geometries, embedded into an elastic or poroelastic medium. Such a fully nonlinear inversion has led to interesting results in volcanic environments and in the study of water tables variation due to its exploitation. This methodology can be used to invert geodetic remote sensing data or terrestrial data alone, or in combination.

Keywords:

inversion technique; 3D models; free-geometry; surface deformation; earth observation data; terrestrial geodetic data; hazards

1. Introduction

Over the past few decades, an explosion in the quantity and quality of geodetic data has been increasingly applied to geological and anthropogenic hazards [1,2,3,4,5,6,7]. It includes ground-, air- and space-based observations that span a range of spatial and temporal resolutions. The development of space geodetic techniques has been particularly important to the measurement of surface displacements, resulting in a sharp increase in the number of active areas where deformation has been detected [5,6,7,8,9]. These new sensitivities and capabilities have produced the need for new approaches for analyzing, modeling and interpreting geodetic data [4,5].

Simultaneously, the rapid growth of the global population concentrated in large cities has led to significant increase in risks when such highly populated areas are also subject to natural (e.g., flood, volcanic and seismic activity) and anthropogenic hazards [4,5,6,7,10,11,12,13,14]. Improving our understanding of the processes associated with these geological and anthropogenic hazards is crucial for developing systems that monitor and evaluate them, as well as to aid decision makers during a crisis.

In the case of volcanic hazards, detecting the accumulation and ascent of magma can provide advance warning of future eruptive activity [15]. Such magmatism is often associated with ground deformation and gravity change, making geodesy a basic tool for the monitoring of volcanic hazards [2,4]. Volcano geodetic observations also provide valuable data for exploring the geometry and volume of magma plumbing systems, which offer a critical context for interpreting the mechanisms and characteristics of unrest and eruption [16,17].

Additionally, land subsidence, ranging from local collapse to the broad regional lowering of Earth’s surface, represents the main geomechanical effect related to the removal of subsurface support [5]. Subsidence can occur as a result of (i) natural factors (e.g., tectonic activity, self-consolidation of recent sedimentary deposits, oxidation and shrinkage of organic soils) [18,19], and (ii) anthropogenic processes (e.g., groundwater pumping [20,21,22,23], urban development [24], hydrocarbon or mining exploitation [25,26]). Land subsidence associated with overexploitation of aquifers is a hazard that commonly affects large areas worldwide.

In both cases, volcanic activity and subsidence associated with aquifer exploitation, the link between data and interpretation is given by direct models and inversion techniques, such that both are critical components to any geodetic monitoring study.

Typically, modeling of deformation and gravity change follows an inverse approach using analytical solutions for volume change or shear faulting within an elastic half space-computationally efficient methods that provide insight into the depth, geometry, and strength of subsurface magma bodies and slipping faults [27,28]. Normally, regular geometries (point sources, disks, prolate or oblate spheroids, etc.) are assumed at the initial stages [27] and the resulting inversion is carried out in a linear context.

In a volcanic context, surface displacements are inverted to infer valuable constraints on the active magmatic sources [2,17,29,30]. Surface deformation is a direct consequence of the dynamics of volcanic plumbing systems, and reflect the shape of magma intrusions, the volume of intruding/arising magma, and the emplacement mechanisms.

Advances in computing power have facilitated rapid development of numerical methods that are capable of both incorporating material and structural heterogeneity [31,32] and investigating the characteristics of the subsurface [33]. They can even be used in inverse approaches [34], blurring the lines between the traditional “inverse” and “numerical” categorizations [5].

Models that combine a suite of geodetic observations can provide exceptional insight into the properties of magmatic systems. Joint assessment of, for example, tilt, GNSS, and differential interferometric synthetic aperture radar (DInSAR) data, leverages the strengths of different data types (in this case, temporal resolution, 3D displacements, and spatial resolution, respectively) to better image magma storage and transport pathways [4].

Some of the most compelling combinations of various types of geodetic data involve joint modeling of deformation and gravity change [35,36]. Typically, gravity and vertical deformation are expected to correlate (accumulating magma results in both mass and volume increases, and vice versa), but a growing body of literature documents instances of gravity change in the absence of deformation [36,37]. Such results may be caused by magma filling or draining from void space [38], magma compressibility [39,40], source geometry and depth [40], or topographic effects [41]. In addition, providing a window into the characteristics of magma plumbing systems, joint modeling of deformation and gravity change can constrain the physical properties of the magma.

Similar (and sometime the same) models can be used [5] to study other kind of geological and anthropogenic hazards, such as land subsidence related to groundwater pumping. This represents a hazard commonly affecting large areas worldwide, often associated with the increasing demand upon groundwater resources due to expanding metropolitan and agricultural areas in semiarid and arid regions [21]. The surface ground deformation thus constitutes a signature of the processes in the reservoir and can provide information about those subsurface processes [42]. Many recent studies have focused on this topic.

The modelling of surface deformation patterns can provide significant insights into the temporal changes of pore pressure as well as the 3D geometry of a reservoir in response to its exploitation over the time [43,44]. A number of different techniques have been developed in recent decades to estimate the surface deformation pattern related to volume changes in elastic and poroelastic media [23,44,45,46,47,48,49,50]. Inverse modeling is required to achieve success in such an endeavor [1,41]. A proper understanding of the subsidence mechanism is essential to calibrate protocols and best practices for monitoring natural and anthropogenic phenomena, with the aim to reduce vulnerability and risk for infrastructures, economies, natural environments and human life.

Analytical models are used to estimate the amplitude and pattern of surface deformation based on assumptions about the media and perturbation source (e.g., using elastic or poroelastic theory) [17,50,51]. They provide a relatively simple method to model surface deformation for reservoirs of any geometric shape. Furthermore, given that these techniques assume that most of the surface deformation is explained by the poroelastic expansion or contraction of the reservoir, less in situ geological data is required than that needed for numerical models [50].

Usually, analytical magma source models represent mathematical abstractions based on strong assumptions (e.g., a priori definition of the source shape, use of geometrically simple geometries, and selection of medium characteristics) that make the set of differential equations describing the problem tractable. Sometimes the analytical approach must be limited to a very simplified and rough approximation of the real source. Nevertheless, analytical models may be useful in certain applications where priority is given to the fast computation of inverse solutions.

On the other hand, numerical Finite Element (FE) modeling allows realistic features such as topography and crustal heterogeneities to be included. However, FE modeling requires large computational resources and it is still very time-expensive to solve the inverse problem with FE models for real time monitoring purposes, i.e., it can take several hours of computation to invert data in order to estimate the source parameters [31,52].

In an operative real time monitoring system, minimizing the delay in the estimation of location and characteristics of the shallow magma source and their evolution in time is the desired goal. Minimum delay in inferring the possible location and timing of an impending eruption would provide critical time for decision makers. Indeed, the fast response of a magma source model can be crucial in early warning for mitigating the possible dramatic consequences of an eruption near populated zones. In this regard, the use of an analytical model for data inversion could speed up the source estimation. Nevertheless, it is still necessary to impose assumptions on the a priori source geometry. Unfortunately, this selection often is highly unconstrained due to the complexity and lack of knowledge of volcano plumbing systems.

Camacho et al. [1] developed an original methodology aimed at the determination of the 3D geometry and location of the causative bodies by inverting ground deformations and gravity changes due to pressure and/or mass anomalies embedded into an elastic medium. Such a fully nonlinear inversion has led to interesting results in volcanic environments, where ground deformations are related to over-pressured magmatic bodies [17,53,54] or water extraction [5]. This methodology can be used to invert geodetic remote sensing data or terrestrial data alone, or in combination.

In this work, we present a review of this method for the simultaneous, nonlinear inversion of gravity changes and surface deformation using bodies with a free geometry. Assuming simple homogenous elastic or poroelastic conditions, the method determines a general geometrical configuration of pressure and density sources. These conditions do not completely represent the complex geological conditions of active areas but are do obtain good results as described in the literature. These sources are described as an aggregate of pressure and density point sources, fitting the entire data set (given some regularity conditions). The approach works in a step-by-step growth process that allows us to build very general geometrical configurations.

2. Inversion Approach for Deformation and Gravity Changes

2.1. System of Nonlinear Equations

Camacho et al. [1] and Fernández et al. [5] consider a set of observation points at the surface (e.g., a geodetic network composed of several stations) (x_l, y_l, z_l), l = 1,…, n, where gravity changes dg_i, i = 1, …, n_g, and position changes (dx_j, dy_j, dz_j), j = 1,…, n_p, have been observed, nearly simultaneously, but not necessarily in the same locations.

Let us denote by d the n_g+n_p observation vector of components dg_i, dx_j, dy_j, dz_j. They try to model this vector by using m buried point sources, located in (X_k, Y_k, Z_k), k = 1, …, m, and characterized by corresponding values for volume v_k, (positive or negative) incremental pressure p_k, and (positive or negative) incremental density ρ_k changes.

They fit the data d as

d = d^c + ε

(1)

where d^c represents the n_g+n_p vector of modeled values

d g_{i}^{c}

d x_{j}^{c}

d y_{j}^{c}

, and

d z_{j}^{c}

for gravity and position changes, and ε is the n_g+n_p vector for residual values coming from inaccuracies in the observation process and also from insufficient model fit [1].

They simultaneously calculate the modeled changes,

d g_{i}^{c}

d x_{j}^{c}

d y_{j}^{c}

, and

d z_{j}^{c},

by means of simple formulas for vertical attraction and deformation effects that are due to additional mass and pressure of several point sources (X_k, Y_k, Z_k; v_k, ρ_k, p_k) within a homogeneous elastic halfspace characterized by a shear modulus μ (given in stress units of pascals, Pa) and a Poisson’s ratio of ν ≅ 0.25.

The surface deformation dx, dy, dz, due to a buried point source is considered mainly as composed by (1) the effects that are due to the incremental pressure p_k and expansion radius within the elastic half-space, described by the Mogi model [31], and (2) the effects produced by the loading of the additional point mass v_k × ρ_k within the elastic half-space, described by the Boussinesq model [55,56].

The surface gravity changes dg can be decomposed as follows [1]:

Free-air effects, corresponding to the relocation of the benchmarks, due to elevation changes according a free-air vertical gravity gradient (about −290 μGal/m for Campi Flegrei); this effect can be included in the model fit using modeled or observed elevation changes;
Newtonian effects due to density changes within the original boundaries of the deep bodies;
Newtonian effects due to mass relocation or change of volume;
Effects due to mass uplift in the surface corresponding to elevation changes. These effects can be obtained using another vertical gravity gradient, depending on the regional terrain density (similar to the Bouguer correction [57,58,59,60]). Effects 1 and 4 depend on the surface elevation changes and can be combined by using a combined gravity gradient F (about −210 μGal/m);
Water table effects which correspond to very local and shallow perturbations;
Topographical effects used for the approximate approach of assuming different source depths at each point [61,62].

Furthermore, they assume possible common and constant terms dx₀, dy₀, dz₀, for the position changes (for instance, an unknown change of the leveling origin) and a possible constant term, dg₀, to represent the simplest form of a long-wave component, a global offset value, or a common uncorrected groundwater effect. Then, with the conditions described previously, they adopt the following expressions for the modeled changes [1]:

d z_{j}^{c} = - \sum_{k = 1}^{m} \frac{g_{0}}{4 π μ d_{j k}} [2 \cdot (1 - ν) + \frac{{(z_{j} - Z_{k})}^{2}}{d_{j k}^{2}}] v_{k} ρ_{k} + \sum_{k = 1}^{m} \frac{1 - ν}{4 π μ} \frac{z_{j} - Z_{k}}{d_{j k}^{3}} 3 v_{k} p_{k} + d z_{0}, j = 1, \dots, n

(2)

d x_{j}^{c} = - \sum_{k = 1}^{m} \frac{g_{0} (x_{j} - X_{k})}{4 π μ d_{j k}} [\frac{(z_{j} - Z_{k})}{d_{j k}^{2}} - \frac{1 - 2 ν}{d_{j k}^{} + z_{j} - Z_{k}}] v_{k} ρ_{k} + \sum_{k = 1}^{m} \frac{1 - ν}{4 π μ} \frac{x_{j} - X_{k}}{d_{j k}^{3}} 3 v_{k} p_{k} + d x_{0}, j = 1, \dots, n

(3)

d y_{j}^{c} = - \sum_{k = 1}^{m} \frac{g_{0} (y_{j} - Y_{k})}{4 π μ d_{j k}} [\frac{(z_{j} - Z_{k})}{d_{j k}^{2}} - \frac{1 - 2 ν}{d_{j k}^{} + z_{j} - Z_{k}}] v_{k} ρ_{k} + \sum_{k = 1}^{m} \frac{1 - ν}{4 π μ} \frac{y_{j} - Y_{k}}{d_{j k}^{3}} 3 v_{k} p_{k} + d y_{0}, j = 1, \dots, n

(4)

d g_{i}^{c} = \sum_{k = 1}^{m} \{G \frac{z_{i} - Z_{k}}{d_{i k}^{3}} - \frac{g_{0}}{4 π μ} \frac{F}{d_{i k}} [2 \cdot (1 - ν) + \frac{{(z_{j} - Z_{k})}^{2}}{d_{i k}^{2}}]\} v_{k} ρ_{k} + F \sum_{k = 1}^{m} \frac{1 - ν}{4 π μ} \frac{z_{i} - Z_{k}}{d_{i k}^{3}} 3 v_{k} p_{k} + d g_{0} + F \cdot d z_{0}, i = 1, \dots, n_{g};

(5)

where

r_{i k} = {({(x_{i} - X_{k})}^{2} + {(y_{i} - Y_{k})}^{2})}^{\frac{1}{2}}, d_{i k} = {({(x_{i} - X_{k})}^{2} + {(y_{i} - Y_{k})}^{2} + {(z_{i} - Z_{k})}^{2})}^{\frac{1}{2}},

G is the gravitation constant, ν, μ are elastic parameters (Poisson ratio and shear modulus), and g₀ is a mean surface gravity value (g₀ ≈ 9.8 Gal).

Equations (1)–(5) represent linear equations for ρ_k and p_k, but a nonlinear relationship with respect to the geometrical parameters (point sources “filled” with non-zero anomalous density or/and non-zero anomalous pressure). We observe that they calculate the free-air and “Bouguer” effects of the elevation changes as proportional to the modeled values

d z_{j}^{c}

, instead of proportional to the observed ones. This allows us to include gravity and elevation changes corresponding to different stations and gives better final results.

In the case of DInSAR data, deformation ds on the surface points is measured along the radar line of sight (LOS) corresponding to the direction of the satellite antenna. Then, the modeled changes

d s_{i}^{c}

for the n data points can be written as

d s_{i}^{c} = d z_{i}^{c} c o s β + d x_{i}^{c} s i n β c o s α - d y_{i}^{c} s i n β s i n α, i = 1, \dots, n

(6)

where α, β are the direction angles (azimuth and incidence) for the antenna pointing direction.

The described system is valid for a volcanic loading problem [1], but if we want to model displacements produced by water extraction in an aquifer [5] we have to move from an elastic to a poroelastic medium. Then, considering the linear theory of poroelasticity [63], the horizontal and vertical components (du, dv, dw) of the movements at a point (X, Y, Z) of the free surface, due to a differential nucleus, located at (x, y, z) and with sides ∆x, ∆y, ∆z, corresponding to the reservoir with local overpressure Δp, are [1,45,46]:

(\begin{matrix} d u \\ d v \\ d w \end{matrix}) = ∆ p \frac{1 - ν}{π} c_{m} (\begin{matrix} X - x \\ Y - y \\ Z - z \end{matrix}) \frac{∆ x ∆ y ∆ z}{{({(X - x)}^{2} + {(Y - y)}^{2} + {(Z - z)}^{2})}^{\frac{3}{2}}}

(7)

where ν denotes Poisson’s ratio (normally estimated at ~0.25) and c_m is the uniaxial compaction coefficient. Assuming that displacements at the surface are almost directly proportional to the thickness Δz of the reservoir, the volume integrations for a parallelepiped cell of sides Δx, Δy, Δz and overpressure Δp in equation (2) can be simplified to integration in the horizontal plane only given rise to [5,46]:

(\begin{matrix} d u \\ d v \\ d w \end{matrix}) = ∆ p \frac{1 - ν}{π} c_{m} I ∆ z

(8)

where:

I = I_{i} (\begin{matrix} X - (x + \frac{∆ x}{2}) \\ Y - (y + \frac{∆ y}{2}) \\ Z - z \end{matrix}) - I_{i} (\begin{matrix} X - (x + \frac{∆ x}{2}) \\ Y - (y - \frac{∆ y}{2}) \\ Z - z \end{matrix}) - I_{i} (\begin{matrix} X - (x - \frac{∆ x}{2}) \\ Y - (y + \frac{∆ y}{2}) \\ Z - z \end{matrix}) + I_{i} (\begin{matrix} X - (x - \frac{∆ x}{2}) \\ Y - (y - \frac{∆ y}{2}) \\ Z - z \end{matrix})

(9)

The integrals I_i for the displacements in the i-direction are:

I_{z} (p, q, r) = \frac{1}{2} \frac{p}{|p|} \frac{q}{|q|} \{a r c \sin \frac{p^{2} q^{2} - r^{2} (p^{2} + q^{2} + r^{2})}{(p^{2} + q^{2}) (q^{2} + r^{2})} + \frac{π}{2}\}

(10)

I_{x} (p, q, r) = a r c \sinh \frac{p}{\sqrt{q^{2} + r^{2}}}

(11)

I_{y} (p, q, r) = a r c \sinh \frac{q}{\sqrt{p^{2} + r^{2}}}

(12)

This formulation provides the direct calculation of the surface effect of a single parallelepiped cell. The total effect of an anomalous structure described as aggregation of m small parallelepiped cells is obtained, according [45], as the summation of the partial effects. This direct formulation can be used to carry out the inverse approach ito determine the 3D pressure source structure responsible of the observed surface deformations for the aquifer case.

2.2. Misfit Conditions

Camacho et al. [1] assume a Gaussian distribution of the observation uncertainty, described by a covariance matrix Q_D for the gravity and position data dg_i and dx_j, dy_j, dz_j (or ds_j), then the minimization condition for residuals ε,

ε^{T} Q_{D}^{- 1} ε = m i n

, leads to the maximum likelihood solution. It is also assumed that the data are not correlated, and Q_D is a diagonal matrix of estimated variances corresponding to observation of gravity and position changes. Gravity and elevation observations are independent and, then, Q_D can be decomposed into two covariance matrices Q_G (n_g × n_g) for gravity and Q_P (n_p × n_p) for position changes. Then the minimization condition for observation residuals can be written as

ε_{G}^{T} Q_{G}^{- 1} ε_{G} + θ ε_{P}^{T} Q_{P}^{- 1} ε_{P} = m i n

(13)

where

ε_{G}

and

ε_{P}

are n_g and n_p vectors for gravity and elevation residuals and θ is a weighting factor for the balance between gravity and deformation fits. This factor is introduced to allow a more versatile fit, making it possible to give optional priority to one kind of observable data over the other. Values for θ can be adopted according to the accuracies of the data sets and considering the approximate equivalence of 1 cm ≈3 μGal.

As is typical for the inversion of geophysical data, problems of singularity and instability for the solution can arise because of inadequate data coverage (e.g., often the number of data points is smaller than the number of unknowns), because of inaccuracy of the data, and because of intrinsic ambiguity of the problem (i.e., if we assume that positive and negative anomalous density and pressure can be contemporaneously present in the model). A possible approach for solving these problems is to use fuzzy logic [64] based on a general estimation of the properties of the physical environment. A more classical and simple process to avoid the instabilities of the solutions is obtained by means of additional minimization or smoothing conditions for the norm of the solution model as

m^{T} Q_{M}^{- 1} m = m i n

(14)

where the model vector m is constituted by the values ρ_k, p_k, k = 1, …, m, for the cells of the model and Q_M is a suitable covariance matrix corresponding to the physical configuration of cells and benchmarks. This matrix provides a balanced model, avoiding very shallow solutions. We use a normalizing diagonal matrix Q_M with elements q_k, k = 1, …, m, given for volumes ν_k and distances d_jk as

q_{k} = \frac{ν_{K}}{n} \sum_{j = 1}^{n} \frac{|z_{j} - Z_{k}|}{d_{k j}^{3}}

(15)

Condition (14) is a stabilizing term [65,66] to control the incremental mass and the pressure for the bodies (weighting according to matrix Q_M). It prevents very large fictitious values of mass and/or pressure resulting from a rather poorly determined model (for instance, one that is due to the coupling of some positive and negative sources, aligned stations, peripheral sources).

A mixed minimization equation,

S (m) = ε^{T} Q_{D}^{- 1} ε + λ m^{T} Q_{M}^{- 1} m = m i n

(16)

is finally adopted as the constraining Equation (1) for residuals and for model magnitude. λ is a factor for selected balance between fitness and smoothness of the model. Low λ values produce a very good data fit (even with noisy inversion) but cam produce an irregular and/or extended model. Conversely, high λ values produce concentrated and smooth models but often with a poor data fit, potentially with autocorrelated residuals. The optimal choice is determined by an autocorrelation analysis of the residual values, as that value producing a null (planar) autocorrelation distribution. See [67] for more details.

2.3. Exploration Approach for Solving the System

The model system, (1)–(5) or (8)–(12), must be satisfied within the minimization constraint (16). The system constitutes a nonlinear problem of optimization with respect to the geometrical properties. Iterative or explorative methods are mostly used for this purpose [68]. Iterative methods (gradient methods) are supported by a long tradition, but they generally require starting with a good initial solution. The more versatile explorative methods can be advantageously applied when the model space to be explored is reasonably sized (small number and narrow range of parameters). In this case, random explorative processes (e.g., genetic algorithm or simulated annealing) are frequently adopted. However, taking into account the very large number of degrees of freedom required to describe the density and pressure models as aggregations of thousands of small cells filled with anomalous values, a general exploratory inversion approach applied simultaneously for all the cells, would be ineffective. An alternative approach is to build the anomalous 3-D structures by means of a growth process and apply an exploratory approach for each step of the growth process.

We carry out a step-by-step process of growth of the 3-D models, using an exploratory technique to find a new cell to be filled with density and/or pressure and aggregated to the models. Therefore, for the kth step of the growth process, k cells have been filled with the prescribed anomalous values for pressure and density, giving rise to modeled values d^c. For the new (k + 1)th step, we try to find and fill a new cell to fit the system:

d = f d^c + ε

(17)

ε^{T} Q_{D}^{- 1} ε + λ f^{2} m^{T} Q_{M}^{- 1} m = m i n

(18)

where f > 1 is a scale factor to allow for a fit between the anomaly of the provisional (not totally developed) model and the observed anomaly (gravity and position changes). To solve, we calculate the value

e^{2} = ε^{T} Q_{D}^{- 1} ε + λ f^{2} m^{T} Q_{M}^{- 1} m

for each of the empty cells according to an exploratory technique. In the (k + 1)th step we choose, as an optimal cell to be filled, the jth cell, one giving

e_{j}^{2} = m i n

(19)

The process continues until we reach a scale factor value f ≅ 1. In the final result, we arrive, nearly automatically, at 3-D models described as the aggregation of point sources filled with the prescribed anomalous values for density and pressure. These simultaneously fit the observed values for gravity and position changes and maintain a small size for both anomalous models.

Camacho et al. [55,69] and Gottsmann et al. [70] provide some simulation examples showing the suitability of this 3-D inversion approach as applied to gravity data. Results are quite satisfactory. A trend towards producing rounded bodies or other distortions is observed for zones of low-reliability, largely because of the additional condition necessary to solve the underdetermined system. Camacho et al. [1] and Fernández et al. [5] carried out some simulation examples to study the performance of the previously described inversion process, divided into two groups. The first group represents volcanic sources as a sum of pressure and mass sources and the second represents water table variation as pressure change sources. It is observed that the fit between the simulated and modeled structures is quite satisfactory [1,5]. The corresponding inversion structure fits the location, depth, and magnitude, but offers a more rounded geometry because of the imposed smoothing condition. In addition, the sizes and depth are quite representative, although some distortions are observed for very peripheral or very deep areas in the model.

3. Application Results

The interpretation methodology, as described above, has been applied to different and complementary test cases. Some of these are volcanic, focusing on the study of a particular period, a long-term time series, or the real time inversion of deformation data acquired during a volcanic crisis. An additional application is the study of the deformation associated with water extraction to determine the volume of the extracted water.

These test cases show the applicability of our methodology for the interpretation of geodetic remote sensing data alone or combined with other terrestrial ones, using only displacement or deformation data, or in combination with gravity variations.

The obtained results for the interpretation methodology demonstrate their applicability to the study of a variety of natural and anthropogenic hazards with different approaches. We describe the application test cases in the following subsections.

3.1. Modeling of Campi Flegrei Unrest 1992–2000 Using Deformation and Gravity Changes

Camacho et al. [1] applied the methodology to actual data from the caldera of Campi Flegrei (Italy), Figure 1. They used gravity and leveling data for terrestrial gravity benchmarks and LOS deformation observed by DInSAR for ascending and descending passes. All the data correspond to nearly the same time period, 1992–2000. They employed change rates, i.e., mean annual displacement or gravity changes. See [1] for a description of the geological setting, the volcanic history, and previous geophysical models.

The leveling benchmarks employed are only a part of a longer network covering the whole area [71,72,73], referenced to an Istituto Geografico Militare Italiano (IGMI) station located in Naples. The double-run survey meets first-order standards. Uncertainties in elevation differences are less than ±2 mm D^1/2, where D is the distance along the line between benchmarks [70], and a final uncertainty of the elevation value d_z at each benchmark result is typically less than 1 cm [74].

In the time interval detained in this work, the gravity network consisted of 15 stations (Figure 1a), each one coinciding with benchmarks of the leveling network. Gravity-height data have an inverse linear correlation [73,75]. As an example, Figure 1b shows the temporal changes for gravity and elevation at the Serapeo station, since it is one of the oldest stations and is close to the area of maximum vertical movements. Camacho et al. [1] observe that the elevation change data show a generally continuous feature in time, but the gravity data contain a higher noise level, with a larger oscillating feature. Then, instead of the real observed data (the values just observed for 2000 and 1992, independent of the other values), they use some filtered (or interpolated or extrapolated) values. This approach allows them to avoid local or instantaneous perturbations or errors in the data, providing a better definition of the inversion model. In addition, they use a mean deformation rate (mGal or centimeters per year) for each station as input data for the inverse approach. See [1] for additional details on leveling and gravity.

For the same area, Campi Flegrei, Small Baseline Subset (SBAS) results (linear velocity and displacement time series) were obtained from a set of 165 European Remote Sensing Satellites 1 and 2 (ERS-1 and ERS-2), and 62 ENVISAT SAR data acquired between 1992 and 2008 on ascending (track 129, frame 809) and descending (track 36, frame 2781) orbits [76]. Camacho et al. [1] use a temporal subset of this data set overlapped with the common time span of the leveling and microgravity data. In particular, data were limited to the periods 1992–2000 for the descending pass and 1993–2000 for the ascending pass.

Ascending and descending deformation linear rate maps are shown in Figure 2. Both maps show a roughly circular pattern of subsidence (increase in line-of-sight phase change) with maximum average velocities larger than 3 cm/yr. Despite some minor uplift [1], the linear velocity during the whole period 1992–2000 is negative. Linear ground deformation rate and time series evolution have been extensively validated against independent geodetic techniques [77]. The global precisions of the SBAS results were determined to be ±0.5 cm and ±0.1 cm/yr for the linear ground deformation rate and time series displacement evolution, respectively [78]. See Camacho et al. [1] for more details in the DInSAR data sets.

In summary, the input data sets used [1] are: (1) The terrestrial data set in terms of linear deformation velocity, including 15 benchmarks with coordinates (x, y, z), gravity change rates dg, and elevation change rates dz for the period 1992–2000; and (2) the total DInSAR data (approximately 15,000 points). For a faster process they selected pixels separated by a distance larger than, for instance, 300 m. It produces two subsets, each with approximately 1330 data points. They apply the inversion methodology to the combined subsampled data set (about 2670 data points) with both ascending and descending data and the terrestrial gravity and elevation, using all the data simultaneously.

A suitable relative weighting for the different types of data is necessary to carry out a simultaneous fit [1]. Leveling data have a good confidence level (better than 1 cm), and their temporal evolution is quite stable (e.g., Figure 1). Gravity data also have a good confidence level (better than 10 mGal), but their temporal evolution is noisier (Figure 1). For the DInSAR data, we have a rather good accuracy level (assumed to be ±0.5 cm and ±0.1 cm/yr). In accordance with these considerations, we assume a similar relative weighting factor for our various types of data (gravity, leveling and DInSAR). As a result, the modeling process is going to be conditioned by these weighting criteria.

The assumed values for the elastic parameters are Poisson ratio = 0.25 and shear modulus = 10 GPa [79]. The assumed gravity gradient used is −220 mGal/m. It corresponds to a free-air value of −290 mGal/m [72] and a value 70 mGal/m for the Bouguer correction for Campi Flegrei, resulting from the analyses carried out by Berrino et al. [80].

To obtain the inverse models they consider a 3-D partition of the local subsurface volume into 14,500 elemental point sources (located at the center of the cells with sides of about 120 m). Using this distribution of “empty” elements, the inversion approach produced the models for pressure and density as defined by the aggregation of filled point sources according to the prescribed extreme values for anomalous density and pressure. After some trials, they selected the extreme values of ±20 kg/m³ and ±10 MPa as a suitable contrast to obtain extended anomalous bodies [1]. The standard deviation of the inversion residuals is 0.2 cm/yr for the elevation changes, 0.6 mGal/yr for the gravity changes, and about 1.1 cm/yr for the LOS values. Figure 3 shows the 3D source model obtained, and Figure 4 shows some additional views of the inversion residuals for the case of the ascending LOS component.

Discussion

The first main result is that the process does not produce any significant volume for anomalous positive or negative density. The resulting anomalous model for density changes consists of only a very few filled cells located close to the surface below some benchmarks, which correspond to very local effects. The anomalous 3-D depressurized model (Figure 3) [1] shows a clear structure located below Pozzuoli at a mean depth of about 1500 m bsl. The vertical profiles in Figure 3 suggest a “partially filled parabolic glass” shape along a WNW–ESE course, with a less-developed similar shape along the orthogonal course. The picture also suggests some small tracks from the main body toward the surface.

In this application case [1], the anomalous mass model obtained from a joint inversion of gravity and leveling changes is nearly empty. There is no anomalous mass below the central area. Only some very shallow (depth about 200 m) and small anomalous (positive and negative) masses are adjusted close to some benchmarks, corresponding to local effects. It suggests that mass changes are not detected below the main deformation area. The gravity changes observed in the benchmarks of the main deformation area can be fully interpreted as produced by free air and Bouguer effects of the uplift (according the assumed gravity gradient). They [1] conclude that the phenomena below the central area do not involve input or output of significant mass (magma or fluids).

The absence of a significant mass change would imply no temporal variability of the fluids in the source region [1], which might be consistent with an exponential decay of fluid pore pressures during the analyzed slowdown subsidence phase (post-1987), following the rapid uplift phases (1982–1984) during which high pore-pressure conditions were established. A physical mechanism that could reconcile these observations could be a long-lasting transient pulse of pore-pressure diffusion confined to the caldera filled material, which does not necessarily imply a mass transfer process (fluids in and out) in the hydrothermal system [81]. However, the absence of mass change does not preclude high pore-pressure conditions from producing fluid flow at the more permeable, shallow, hydrothermal system levels, as suggested by the geochemical variation of volatile discharge at Solfatara [82]. Indeed, local fluid flow could explain the few very shallow mass change sources required to fit some benchmarks.

All deformation data (terrestrial, ascending, and descending DInSAR) provide a similar low-pressure structure located at a shallow depth (about 1500 m) below the survey area [1]. The shape of that structure nearly fills the bottom of a parabolic cup with the same additional filling in the walls (Figure 3). Therefore, they did not detect signs of the deep magmatic source for this period. Some main WNW–ESE elongation is detected at about N100°E. The parabolic cup is superimposed with depth on the shape of some existing structure, namely the inner caldera edge, as depicted by a 2.5-D inversion of on-land and offshore gravity data [83]. Moreover, its bottom corresponds to the depth of the bottom of the inner caldera (between 2 and 2.5 km). The modeled structure fits well the polygenic body inside the caldera. We interpret this as corresponding to the dynamics of the shallow (depth 1–2 km) hydrothermal system, as previously suggested by several authors [74,84,85]. The hydrothermal system is confined to the filling caldera materials and is limited by the inner caldera structure.

3.2. Modeling of Campi Flegrei Unrest 1993–2013 Using Only Displacement Data

DInSAR is used extensively today for mapping ground deformation with high spatial resolution and sub-centimeter precision over large areas and is a suitable tool for ground deformation monitoring of active volcanic areas [4,53]. Spatial resolution of modern synthetic aperture radar (SAR) sensors ranges from 1 to 20 m over areas from 10 × 10 km² to 250 × 250 km² [53]. For modern satellite constellations the repeat cycle ranges from 1 day to a few weeks, but the typical repeat cycle for a single satellite mission ranges from 12 to 41 days. Repeatedly acquired SAR data from a single sensor are used to obtain line-of-sight (LOS) time series of ground deformation by applying small baseline subset (SBAS) [86,87,88], persistent scatterer [89] methods, or a combination [90]. The results of these techniques, however, are limited to the time period of the particular data set and do not distinguish horizontal and vertical motion.

The multidimensional SBAS (MSBAS) technique [53,91] combines multiple DInSAR data sets into a single solution with improved characteristics: Lower noise, denser temporal resolution than earlier DInSAR time series, and continuous temporal coverage. The MSBAS methodology is a natural extension of the original SBAS method that addresses the data redundancy and multidimensionality of the problem by decomposing LOS DInSAR measurements into vertical and horizontal (mainly east-west) time series of surface deformation using ascending and descending DInSAR data. MSBAS has been applied to the mapping of anthropogenic [92,93] and natural [91] ground deformation, successfully producing 2-D time series with dense temporal resolution and high precision.

Samsonov et al. [53] apply the MSBAS DInSAR processing technique to 20 years of ERS-1/2, Envisat, and RADARSAT-2 data from Campi Flegrei. The resulting spatiotemporal signals were determined with a duration and resolution that had not been achieved before in this or any other volcanic region. The used radar data are listed in Table 1. The theoretical derivation of the MSBAS technique is described in detail in [53,91].

MSBAS processing starts after both ascending and descending images are acquired [53], so the time series begin on the first date of the ascending ERS data (10 January 1993). The resulting 1271 highly coherent interferograms were used in the MSBAS processing, producing vertical and east-west time series spanning 20 years (1993–2013) and 385 time steps. The average error in the east-west time series is approximately 0.07 cm, while the error in the vertical time series is approximately 0.09 cm, although this varies spatially and can be as much as twice that in regions of maximum deformation (supporting information in [53]). In regions such as Campi Flegrei, where there is a significant trend, this estimate provides the worst case scenario, and the actual error is likely smaller. See [53] for more details on the DInSAR processing.

Figure 5 presents the results of the MSBAS processing [53]. Figure 5a shows the vertical change in surface height between the initial and final time steps. The associated deformation time series is shown in Figure 5c for one location and presents a more complicated picture than the net subsidence of Figure 5a,b (deformation time series for selected locations are presented in the supporting information by [53]).

The precise deformation time series captures more than 3 cm/yr of maximum subsidence (green dot) between 1993 and 2000, centered on the caldera. Subsidence continued at a slower rate, interspersed with short periods of uplift, until 2005. Almost continuous uplift began in 2005 and accelerated to a rate of approximately 2.5 cm/yr between 2008 and 2013. Deformation is ongoing at a rate of 5 cm/yr (2011–2013). The large number of time steps and precise measurements are evident in both the vertical and east-west time series (Figure 5c). Additional times series can be seen in [53].

Using this new data set, Samsonov et al. [53] determine the potential sources below the caldera using the described methodology [1]. They apply the inversion method to specific time periods related to the inflation and deflation at Campi Flegrei observed in Figure 5c, 1993–2013. Table 2 shows those time periods and the resulting parameters for the associated inflation or deflation sources. The various periods of subsidence and uplift (deflation and inflation) can be modeled as separate sources. The results are summarized in Table 2 and in the next Figures.

We show the results from the nonlinear inversion method for different time periods (1993–1999, 1999–2000, 2000–2005, 2005–2007 and 2007–2013) and the actual and modeled surface deformation rates (cm/yr) for both vertical and east-west motion for each of the considered time intervals together with the residuals between the observed and modeled displacements. See Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15. The models resolve the actual displacements quite well, with the exception of coherent residuals in the Solfatara region that are likely associated with the shallow source under the fumarole field [94] and not modeled in that work.

Discussion

The computed deformation by Samsonov et al. [53] show consistent caldera subsidence for the period 1993–1999. However, deformation changed from subsidence in 2001 to primarily uplift in 2005, in conjunction with changes in geochemistry indicators such as CO2/H2O ratios [95], indicating an increase in relative enrichment of a magmatic CO2 component in sampled fumarolic gases. The rate of uplift continued to increase through the present, reaching more than 5 cm/yr.

Comparison of the deflation and inflation sources for the various activity periods shows that they are remarkably similar in location, shape, and volume (Table 2 and Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14 and Figure 15). The inflation source tends to be deeper, nearer the bottom of the marine precaldera and postcaldera deposits and limited in depth by the transition to the gas-bearing deposits layer [53]. The deflation source is generally shallower, near the top of the marine precaldera and postcaldera deposits. This behavior suggests that inflation-deflation cycles are generated by sources fed from below and deeper than the 3 km detected here. Migration of fluids seems to be focused at the top of the gas-bearing formation, overcoming lithostatic pressures and effectively transferring stress to the surrounding rock during inflation periods. This is followed by gradual relaxation of the overpressure (e.g., pore fluid pressure diffusion) at shallower depths [1]. While others have modeled deformation at Campi Flegrei with two simple sources at 1700 m and 3600 m in depth [94], here the bottom of the inflation source extends below 3000 m. In this case, one extended source model with unconstrained geometry explains the deformation [53]. Samsonov et al. [53] showed that the deformation driving mechanism is strongly controlled by the caldera structure and stratigraphy and must involve large extended sources in a layered hydrothermal system. Due to the location of Campi Flegrei and its proximity to the city of Naples, more active monitoring for signs of possible imminent eruption is warranted.

Apart from the geological and geophysical conclusions, their study demonstrates that a temporal resolution on the order of one or two weeks can be achieved for MSBAS time series, nearing that of continuous GPS (cGPS) daily time series, by simultaneously processing DInSAR data from different sensors. The accuracy of MSBAS vertical and east-west components is very high, potentially higher than cGPS. At the same time, the spatial resolution of MSBAS results is superior to cGPS. Therefore, the MSBAS approach can be used for monitoring volcanic activities in near real time, particularly in remote areas where other observation techniques are unavailable or insufficient.

3.3. Real Time Tracking of Magmatic Intrusions During Volcanic Crisis: 2008 May Etna Eruption

The goal followed by Cannavò et al. [17] is to estimate in real time, for a given volcano showing signs of unrest, the location, size and volume distribution of the deformation source and to track those parameters through time without any a priori assumption on the source geometry. They seek to delimit the area with the highest likelihood of being the location of an opening fissure or vent for the impending eruption. They apply the approach by Camacho et al. [1] to model volcanic sources and track their evolution in time. Their objective is to address the limitations in the inversion process dictated by the running time and the a priori assumption of the shape of the intrusion or magma chamber in order to produce a useful real time that only adds seconds to the time needed to obtain displacement data. This methodology [17] performs a nonlinear inversion of the displacement data, producing extended bodies [1,53] with a free geometry within seconds (less than a minute) in a fully automatic way, i.e., without requiring any user input.

To demonstrate the applicability and results obtained using this methodology in real time tracking of magmatic intrusions using ground deformation data they [17] select the main eruption at Mt. Etna, Sicily, Italy, in 2008 as a test case. The 13 May 2008 flank eruption on Mt. Etna occurred on the upper eastern flank after an intense recharge phase beginning in early 2007, interrupted by seven fire fountain episodes at the New South East Crater (NSEC), the last of which took place on 10 May 2008. The intrusion produced large ground displacements captured by the volcano cGPS network [96].

Due to the magnitude of the event and the quality, quantity and variety of available data, this eruption was a good benchmark to test new volcano modeling techniques. In particular, this study models the ground deformation pattern associated to the co-intrusion dyke propagation in the shallower part of the northern sector of Mt. Etna. The high rate GPS (HRGPS) data for 10 stations (Figure 16) from the Etnean GPS permanent network (Etn@net) is employed to illuminate the dynamics of the magma intrusion [97]. A real time scenario is simulated by computing HRGPS solutions from the raw data and running the inversion code, thus simulating a real processing and inversion stream during the eruption process. The real time analysis starts the day before the eruption when GPS signals showed some transients, though it should be noted that the shallow co-intrusion period lasted from the 8.30 am UTC [96] (when the seismic swarm started) on 13 May 2008 to 04.00 pm UTC the same day (when almost all the seismic activity had ended). The real effusion [93] phase, together with lava fountaining, started at 9.30 am UTC, as revealed by Meteosat Second Generation satellite. The fountaining stopped at around 11.00 am UTC [98]. This implies that the early warning window to localize the impending eruption site was only 1 hour (from 8.30 am to 9.30 am UTC).

Cannavò et al. [17] reproduce the real time scenario of the eruption. Though the test data are processed after data collection (post-processing), real time operation was simulated. Displacements are calculated by a robust linear regression [99] of the GPS solutions in the 30-minutes window. An interval of 30 minutes was chosen as a tradeoff between the minimum response time of the warning system and the minimum level of noise in the estimated displacements (Figure 17). See [17] for more details on the GNSS processing method.

Inversions are performed in a continuous way from the simulated real time displacement data computed every 15 minutes. There is some general uncertainty for any inversion approach with respect to the pressure contrast ΔP of the causative structure. Cannavò et al. [17] use 1 MPa as the value of pressure change for inversion in this test case, considering typical values obtained studying different eruptions for Mt. Etna volcano and from other studies of the considered test case eruption [96]. Authors perform several tests, considering different values of pressure change (0.5 to 3 MPa), that produced very similar results.

3.3.1. Results

Before testing the real time response of the inversion procedure, Cannavó et al. [17] use the inversion methodology to analyze the inflation phase preceding the 13 May 2008 eruption. They model the time evolution of the magmatic source for the period 1 January 2007 to 12 May 2008 by considering 15 day-intervals (Figure 16 and Figure 18). It is determined that, during the months of June and July 2007, the system was subject to recharging episodes characterized by an ascending high-pressure body. As shown in Figure 19, the modeled body ascended from a depth of about 10 km, to reach a level of about 2–3 km depth (b.s.l.). This modeling study allowed estimation of the rate of magma ascent and suggests an approximate 3D geometry. It was assumed that an a priori estimated knowledge and interpretation of mid-term (weeks-to-months) ground deformation signals are routinely obtained in most volcano observatories. Such estimated models account for typical monitoring resources, which include a weekly report of volcanic activity, integrating the available geodetic data on those time scales.

Nine months after the deep magma ascent, a dyke intruded in the North-East sector of the Etna summit area. The simulated real time inversion produced a time series of 3D body sources shown in Figure 20. In particular, they observed that in the evening (20–22 UTC, all the times are in UTC hereinafter) of the day before the eruption, a vertically elongated source from 2 km bsl to 2 km asl (above sea level) identifies a magma batch likely migrating towards the surface (see the Supplementary Videos [17]). Starting from 6.00 am of the day of eruption, a large positive pressure source located at about 2.5 km asl in the NE sector of the summit area triggered the dyke intrusion at 7.00 am towards the shallow complex source in the NW-SE direction (Figure 19, panel a). The latter source reached its maximum expansion in that direction at 11.30 am (Figure 19, panel c) and contracted only after 12:30 pm. This activity exactly follows the sequence of events and models reported in previous studies [96,98]. The estimated pressure sources show, minute by minute, that most of the recorded deformation corresponds to a shallow dilation process located close to the NE sector of the summit area, just where the eruption fracture opened.

Figure 16 shows the traces of the major sources that emerged from the inversion for both the studied time scales (long and short) under consideration together with seismicity during these periods. The uncertainty of the source location is estimated at about 500 m in the horizontal component and 900 m in the vertical component [17]. Thus, by considering the associated errors of the earthquakes (200–300 meters in all the 3D components) recorded by the local permanent seismic network managed by INGV, from Figure 16 they can estimate the earthquake swarm to be co-located with the estimated magma sources during the intrusion phase. It is worth noting that all the identified elongated model sources seem to occur near the periphery of the High Vp Body (HVB) [100,101]. A comparison between the shape and the position of the volumes estimated by Cannavò et al. [17] and the ones modelled by Bonaccorso et al. [98] for the period 1993–2000 suggests that there is a magma reservoir at about 2 km bsl and a vertically elongated source located in the south-western border of the HVB, where magma migrates towards the surface. Figure 20 provides a 3D-view sketch summarizing the results.

The goodness-of-fit for the estimated time varying model is assessed by computing the root mean square error (RMSE) between the measured and the modelled displacements. Discrepancies between observed and modelled values are in the order of 0.2 and 0.4 mm/day for horizontal and vertical components displacements, respectively [17].

Some spurious sources appear scattered around in the domain (see e.g., Figure 19 and the Supplementary Videos by [17]). They primarily are y due to the residuals in the data (non-ideal network topology configuration, noise measurements, and errors in the parameter settings) that the algorithm attempts to account for. A suitable trade-off between good data fit (with possible aggregation of spurious sources) and good model stability (with poor data fit and too simple model geometry) is obtained by analyzing the autocorrelation of the resulting data residuals [1]. Nevertheless, spurious sources are usually easily distinguishable, but sometimes they aggregate, forming a volume that may appear as a real source. In such cases, considering the spatial location, the size, the shape and their temporal evolution, they can usually be identified and discarded as true sources of deformation.

In order to quantify the reliability of the solutions obtained with this method, Cannavò et al. [17] carry out an uncertainty analysis. Figure 21a shows the decrease of the sensitivity to pressure with the depth and with the distance to the stations. In fact, the required pressure (for a 1 km³ body) to produce a particular magnitude of surface deformation at the current network stations (root mean square, rms, 1 mm) increases with depth and distance to the stations. For deep and distant pressure elements in the inverse model, the relative uncertainty increases with a pattern similar to that of Figure 21a. Figure 21b,c show the 3D uncertainty distribution for vertical and horizontal displacements of the isolated pressure sources. It also is valid for some geometrical components or details of extended sources.

These maps show the location values (in km), for a source volume 1 km³ and a pressure value 1 MPa, required to produce a particular surface deformation level (1 mm, rms value). They can be interpreted as showing the pattern of relative uncertainty for the geometrical aspects (depth and horizontal location) for the inverse model. In addition, the displacement uncertainty pattern is heavily dependent on the pattern of the uncertainty for the source intensity (pressure and volume) of the assumed source element.

A more valuable uncertainty map for model geometrical aspects is obtained by assuming equally sensitive source elements (according to the pattern of Figure 21a). Figure 21d shows the sensitive 3D pattern for displacements for a volume of 1 km³, but assuming pressure values increasing with distance according to Figure 21a. These last model uncertainties show a different pattern, essentially due to geometrical aspects. Further details on uncertainty analysis can be found in the Supplementary material by [17].

The uncertainty analysis shows the power resolution of the GPS network for the proposed source estimations. The uncertainty maps show that the sources are well-resolved in the volumes where they have been estimated in Figure 16, Figure 18 and Figure 19, thus making the results more reliable.

To analyze the uncertainty of the results to the GPS network configuration, Cannavò et al. [17] carried out the same sensitivity tests reported in the manuscript by using the latest GPS network on Mt. Etna. Here they estimated the results for the current configuration of the GPS network composed of 40 stations. The results show a rather different pattern, mostly corresponding to the geometrical distribution of stations. For some common areas the 40 stations sensitivity is obviously higher (smaller pressure changes or smaller displacements) than that corresponding to 10 stations. This means that the more stations involved in the inflation/deflation process (i.e., measure displacements), the better the resolution of the model [17].

3.3.2. Discussion

The main goal of this test case is to study the ability of the described inversion approach to instantaneously provide a distribution of pressure sources for real time volcano monitoring purposes. Cannavò et al. [17] show that this inversion methodology overcomes the limitations of fixed-shape analytical models and the time-consuming nature of numerical approaches based on finite elements, providing a tool to infer the timing and location of the potential impending eruption, thus helping to activate warning procedures in good time. It may therefore be a crucial tool for civil protection purposes. For example, during the first hours of the studied 2008 Etna eruption, bad weather conditions prevented researchers from identifying the location of the ongoing eruption. However, the use of a tool like the one proposed here could have assisted researchers during the peak of the crisis by helping define the magma location in real time and infer the area of the opening eruptive fracture.

The hypothesis of homogeneous elastic medium does not appear to affect in any critical way the quality of the test results considered here, as was expected with the time periods used [103,104]. They show good agreement with previous results described in literature, e.g., the shallow NW-SE oriented source described above overlaps the previously modeled dyke [96].

Tracking the evolution of the magmatic source beneath the volcano surface in real time has always been a challenge for researchers working in modern volcanic surveillance. The proposed methodology is able to estimate, in real time, the location, size and shape of the active magma source in volcanic areas and can track those parameters through time, without any a priori assumption on the source geometry. This information can be the basis for forecasting activity in a monitoring and hazard assessment system to assist decision makers during a volcanic crisis. To this end, this test case showed the effectiveness of the proposed methodology for one of the best monitored active volcanoes on Earth. Mt. Etna represents a significant hazard to the local population living in neighboring areas and to buildings and infrastructure. The methodology, though successfully applied here to a Mt. Etna eruption, can be tuned and used in all active volcanic areas. Compared to other approaches in volcano monitoring by means of geodetic data, the advantages of the proposed tool are: (1) The speed of data inversion (few seconds); (2) the absence of assumptions on the shape of the source, (3) obtaining a 3D geometry of the magmatic sources and their time evolution; and (4) the easy implementation of the tool (it runs on a common personal computer without any special requirements).

3.4. Land Subsidence Associated with Overexploitation of Aquifers: Lorca, Spain

The Lorca region, located in the Alto Guadalentín Basin of southeastern Spain (Figure 22), is affected by subsidence rates of up to 10 cm/yr as a direct consequence of long-term aquifer exploitation [5,21,22,23,24,25,26,27,28,29,30,31] (Figure 22). This region is characterized by semi-arid climate conditions, with average precipitation rates of 150 mm/yr and an average annual temperature [105] of ~18 °C. The basin is infilled with Quaternary alluvial fan systems overlapping Tertiary sediments transported by the Guadalentín River along the depression located in the eastern part of the Betic Mountain Range (an ENE-WSW oriented alpine orogenic belt resulting from the Nubia-Iberia ongoing convergence [106,107,108]).

The Guadalentín Basin aquifer is composed of two contiguous sub-basins: Alto and Bajo Guadalentín (Figure 22). From a hydrogeological point of view, the basement beneath the aquifer is composed of several relatively impermeable Paleozoic metamorphic complexes overlain by permeable Miocene conglomerate and/or calcarenite series. The top of the succession comprises Pliocene-Quaternary, low-permeability, compressible conglomerates, sand, silt, and clays [21,109]. The Alto Guadalentín aquifer covers an area of approximately 277 km². Historically, piezometric levels were located closely to the land surface, allowing the development of a number of artesian wells and permanent lagoons [109].

Since the 1960s–1970s, the Guadalentín Basin aquifer has reflected gradually increased overdraft and contamination (e.g., high electrical conductivity, CO2 positive thermal anomaly), and was legally declared provisionally overexploited [110] in 1987. Pumping has occurred in ~1000 wells at rates of 24 (in 1973), 69 (in 1987), and 86 hm³/yr (in 2006) [109,110], which led to a spatially variable continuous piezometric level decrease (at rates within the 0.5–10 m/yr range).

Available piezometric information consists of water-level time series for a few points, from the 1970s to present. A long drought period from 1990 to 1995 (also in 1999–2000 and 2005–2007) reduced natural recharge and increased pumping in the Guadalentín Basin, which led to an increased resources deficit. All this information indicates a long-term trend in the consumption of groundwater resources [5,21].

Interferometric synthetic aperture radar (DInSAR) studies, while detecting the high subsidence rates affecting the Alto Guadalentín Basin, also identified a delayed transient of nonlinear compaction of the Alto Guadalentín aquifer due to the 1990–1995 drought period [21]. This suggested a relationship between local crustal unloading and stress change on active faults bordering the basin [51]. Later work [22] extends those studies using advanced differential DInSAR (A-DInSAR) techniques to process ALOS PALSAR (2007–2010) and COSMO-SkyMed (2011–2012) radar images. The combination of multi-sensor SAR images with different resolutions allowed for a longer monitoring time span of 20 years (1992–2012) over the Alto Guadalentín Basin. Additionally, the satellite measurements study area. Furthermore, the work presented a new soft soil thickness map and collected historical piezometric data, in order to assess aquifer system compressibility and groundwater level changes in the past 50 years. From the analysis of these data with A-DInSAR displacement measurements, the authors conclude that the governing mechanism of the Alto Guadalentín aquifer system is an inelastic, unrecoverable and delayed compaction process between water level depletion and ground surface displacement, related to the presence of very thick (>100 m) unconsolidated sediments (clay and silts).

Despite the aforementioned achievements, the previous studies focusing on the deformation in the area are based on DInSAR analysis using ascending and/or descending acquisitions, without any combination of the datasets to estimate both vertical and horizontal (E-W) components [21,22]. Therefore, only the line-of-sight (LOS) displacement field is known in the Alto Guadalentín area at a regional level and it was assumed to correspond completely to vertical displacement. Although this is a common procedure in subsidence studies using DInSAR measurements [112,113,114,115,116,117,118], the main consequences are (i) the neglecting of possible horizontal displacement components and (ii) the likely overestimation of vertical displacement.

Here, with the decomposition of LOS measurements in the E-W and vertical components over the investigated area, additional constraints are provided for the spatial and temporal evolution of the subsidence process as well as on the main governing mechanisms (e.g., temporal changes of pore pressure, geometry of the reservoir). With this primary aim, Fernández et al. [5] established a GNSS network consisting of 33 stations in 2015, which densely covers the Alto Guadalentín basin (Figure 22). This network was observed in survey mode. They analyze the measurements carried out in November 2015, June-July 2016 and February 2017 [5]. GNSS raw data was processed by adopting standard processing strategies for this type of network and referred to a local reference frame in order to estimate the 3D deformation field (see [5] for additional information). Despite the limited time interval covered by the surveys, they estimated, for the first time, a significant 3D deformation field, which is primarily related to the local exploitation of the aquifer.

SAR data from the Sentinel-1 Copernicus constellation, acquired in ascending and descending orbits for the same time period, also were processed to obtain the respective LOS displacements. Using the GNSS and SAR-based deformation fields, both the vertical and horizontal components of the displacement over the entire area were estimated [5]. Finally, they use these results for comparison between them and interpretation using the forward model and inversion technique previously described.

Fernández et al. [5] carried out three Global Navigation Satellite System (GNSS) campaigns, in November 2015, June-July 2016 and February 2017. The 3D velocity field results are shown in Figure 23 for both the vertical and horizontal components determined by comparing the coordinates obtained for the time spans of the three surveys. See [5] for additional details on the observation and processing methodologies, and on the results.

The maximum vertical subsidence rate (9.0 ± 0.5 cm/yr) is of the same order of magnitude as that previously detected by earlier DInSAR studies [21,22]. The maximum horizontal displacement rate detected is 2.5 ± 0.3 cm/yr (about 28% of the vertical displacement rate), a non-negligible amplitude. In the area showing the highest subsidence rate, again previously detected by DInSAR techniques, a characteristic pattern of horizontal deformation appears (Figure 23). These deformations, as theoretically expected, show a centripetal pattern towards the zone of maximum subsidence, located in the central part of the monitored area.

In the southern area, where there is a relative maximum in the LOS displacement detected by DInSAR, significant horizontal motions are detected (Figure 23) associated with GNSS stations 23 and 28. After comparison with A-DInSAR results and field inspection, we conclude that these are produced by very local movements related to monument instabilities [5].

The Advanced Differential Satellite Synthetic Aperture Radar (A-DInSAR) processing of the Lorca area Sentinel-1A images (see [5] for a detailed description of the advanced processing of the satellite radar images). Both orbits, ascending and descending (tracks 103 and 8 respectively) from Sentinel-1A, are used to decompose the measured LOS movement/mean velocities into horizontal (E-W) and vertical components [119,120,121] over the studied area. The A-DInSAR study covers the same time interval spanned by the GNSS campaigns (November 2015–February 2017). Radar data are processed using the Coherent Pixel Technique (CPT) [5,122]. The total area covered by the GNSS network is approximately 70 km². A-DInSAR was processed for an extended region with a total area of 170 km². In both geometries, ascending and descending, the study area is covered by three bursts of the same swath. A total of 42 Interferometric Wide Swath (IW) SLC images from the Sentinel-1A satellite resulted in 185 interferograms (22 images and 137 interferograms for ascending data, 20 and 48 respectively for descending [5]). The results are shown in Figure 23. Some selected time series are shown in the Supplementary Material by [5] for descending LOS. This is the first DInSAR study of the Alto Guadalentin Basin using two different geometries for the same time period.

Fernández et al. [5] obtain the vertical and E-W components of the displacement from the ascending and descending LOS motions. E-W and vertical components of the displacement fields in the area obtained using ascending and descending LOS results are shown in Figure 24. For this case [5], the decomposition above does not produce particularly good results due to methodological aspects and to the fact that the magnitude of the E-W motion is at sub-centimeter levels in many of the coherent pixels, i.e., the same order of the A-DInSAR uncertainty.

The comparison of GNSS and A-DInSAR results allows to estimate the error of the A-DInSAR processing, relative to GNSS, as ~0.7–1 cm [5].

Discussion and Inversion of the 3D Displacement Field

As previously noted, the establishment and observation (spanning November 2015–February 2017) of a local GNSS network allowed, for the first time, for measurement of the 3D displacement field in the Alto Guadaletin area, associated with exploitation of the local aquifer [5]. Also, for the first time, A-DInSAR results have been obtained using both ascending and descending radar images from the Sentinel-1A Copernicus radar satellite, allowing estimation of both vertical and horizontal (E-W) displacement components, a 2D displacement field, at higher spatial resolution than GNSS [5]. See Figure 23 and Figure 24.

These results [5] highlight how the ad hoc establishment of survey mode GNSS networks improves the spatio-temporal monitoring of the 3D displacement field of areas subjected to extensive groundwater extraction, therefore representing a valuable monitoring technique. Moreover, GNSS observation provides complementary information to A-DInSAR results, allowing for their validation and scaling. In addition, at a local level, it is observed that the GNSS network does not completely cover the current displacement area, in particular along the SW region (Figure 23 and Figure 24), because the network was defined based on displacements obtained prior to 2012 (Figure 22). However, their Sentinel-1 A-DInSAR results [5] show that the deformation has extended in the SW direction, which today is the region of the most significant water extraction [123]. Therefore, the GNSS network needs to be extended over that region with additional GNSS stations.

The results obtained for the studied area highlight that [5]: (i) Simultaneous GNSS and A-DInSAR results are consistent with each other (Figure 24 and Supplementary Information from [5]; (ii) the results obtained for rates and pattern of the displacement are consistent with previous DInSAR results [21,22]; however, (iii), the horizontal displacement rate has a maximum amplitude of 2–3 cm/year (Figure 23) and it is a significant component of the observed deformation field. Therefore, the horizontal displacement cannot be neglected, as it was done in previous studies [21,22].

Since the results by Fernández et al. [5] demonstrate that the horizontal displacements represent a significant component of the deformation field of the studied area, they carry out sensitivity tests by neglecting/including this horizontal motion in order to assess the variability (or bias percentage) on the determination of the aquifer characteristics and their temporal evolution using deformation modeling. To do this, they employ both the GNSS and A-DInSAR results. In addition, taking into account the linear time behavior of the displacement field [4], they also consider displacement rates in their study. Fernández et al. [5] use ten different data sets of surface displacement (denoted as Cases) covering the period November 2015 to February 2017 and carry out the inversion using the previously described forward model and inversion methodology. The cases are the following:

(A) LOS A-DInSAR results obtained for descending orbit images, assuming 100% as vertical displacement;

(B) LOS A-DInSAR results obtained for ascending orbit images, assuming 100% as vertical displacement;

(D) Purely LOS A-DInSAR results obtained for descending orbit images;

(E) Purely LOS A-DInSAR results obtained for ascending orbit images;

(F) Purely LOS A-DInSAR results obtained for ascending and descending orbit images;

(G) Up-Down and E-W A-DInSAR results obtained by combining ascending and descending orbit images results;

(H) 3D displacements determined using the GNSS surveys results;

(I) Purely LOS A-DInSAR results obtained for ascending and descending orbit images together with the 3D displacements determined using the GNSS surveys results;

(J) Up-Down and E-W A-DInSAR results obtained by combining ascending and descending orbit images, together with 3D displacements determined using the GNSS surveys results.

Cases A-C are one-dimensional (1D), D–G are 2D (either indirectly by combining Up-Down and E-W in any of the two measured LOS, or directly supplying the values for both displacement components separately), H is a purely 3D dataset, and I and J are a combination of 2D and 3D data (2D+3D data).

A very important difference between the purely GNSS 3D data and the rest of these cases is the number of displacement rate data points [5]. For the 3D Case (H) data set, there are just 108 displacement rates, while for the remainder we have thousands of measurements (Table 3).

Fernández et al. [5] inverted each case and estimate the volume changes of the water table (volume and geometry) assuming a given pressure change value. Moreover, based on hydrogeological observations, they impose the criteria that sources were shallower than one kilometer. A summary of the results obtained is provided in Table 3 and Figure 25. In these figures, the blue colors indicate negative pressure values, while white colors indicate positive pressure changes cells.

Inversion results include the volume and geometry of the active part of the aquifer which produced the measured displacements. Here this is quantified by the intensity, equal to the product of volume by pressure change, as it is impossible to determine both quantities separately. If we increase pressure, we decrease volume and vice versa. In order to determine a general geometry, they constrained the value of pressure change [1,17]. After a trial analysis, they considered a pressure value of –3 MPa [5], selecting the value that gives us a source geometry most consistent with the characteristics of the Lorca aquifer.

The inversion results obtained from these data sets can be organized into two subsets: (i) Cases A–C (1D, vertical displacement) and (ii) Cases D–J (2D and 2D+3D). Both sets of results are internally consistent, with respective scattering on the order of 9 and 6% respectively (Table 3).

The results of (i) are, on average, ~26% greater in intensity/volume than those of (ii), indicating that using only one component of the displacement field and assuming that the displacements are only vertical significantly overestimates the volume of water extracted during the study period (on the order of tens of hm³). This can have an important effect in predictions of future volume variations and surface displacements.

The results for group (i) show the largest discrepancies between Case C and the other two cases (A and B). For Case C we obtain an intensity value ~20% greater than the other ones. This data set has been obtained from combining LOS results from ascending and descending orbits and they suffer from the small systematic errors in the methodology [5]. This problem is reflected in two aspects, the poorer misfit values and the appearance of additional sources at the edges of the study area that try to adjust these errors (Figure 25). As a result, it was concluded [5] that it is clearly more appropriate, in agreement with previous works [42], to mitigate this problem by using the ascending and/or descending LOS directly in the inversion procedure. As a result, Case C is not considered further. Cases A and B obtain very consistent results (Table 3), with a dispersion of ~4% for group (i) and a mean intensity value of 40 MPa∙km³ [5].

For the group (ii) results, the data set for the vertical and E-W components determined using the combination of ascending and descending LOS (Cases G and J) (Table 3) again provides worst results (higher intensities, more additional sources and poorer misfit values, Figure 25 and Table 3). Case H (only GNSS determined displacements) also provides higher misfits (Table 3) and poorer adjusted geometries (Figure 25h), probably a result of using data with several important limitations [5]. The first one is the reduced number of data points available in relation to the number of unknown parameters (approximately one hundred versus several thousands (Table 3). As previously mentioned, in cases like this, a large number of measurement points are needed to reliably estimate the distribution of reservoir volume changes [5,48]. Also, their limited number produces a poorer spatial distribution than A-DInSAR data and does not cover the entire deformed area (Figure 23). Finally, as described previously, in this reduced data set we have some anomalous results associated to very local effects (Figure 23b and Supplementary Information from [5]). To establish a denser and more extended GNSS network requires additional cost and time for field observation.

Another important result [5] is that significant differences (3–4%, at the level of error or lower, see Table 3) are not observed between using just one LOS (ascending or descending), both LOS (ascending and descending) displacement data sets together, or both combined with the GNSS results. The estimation of source characteristics is very similar for all cases, just slightly changing the misfit of the minimum values for Cases D–F, I.

In summary, contrary to previous studies in the Lorca area, Fernández et al. [5] conclude that the measurement and use of the horizontal and vertical displacements at the surface is important for the prediction of future volume variations and surface displacements. These differences can have important effects on the design of monitoring systems, help in the decision-making process related to the sustainable management of the aquifer resources, and improve the assessment of potential hazards related to the aquifer exploitation. These main conclusions did not include any local assumptions which could condition our interpretation methodology [5], but they have general applicability worldwide. The Lorca case can be considered an extreme case, taking into account that it has significant E-W horizontal deformation, but only in geographically limited areas of maximum deformation. In other regions, where significant horizontal E-W deformation may be more scattered and cover more extended areas [5] the effect of considering vertical deformation alone could be even more dramatic.

In summary, Fernández et al [5] show, using inversion results from different remote sensing data sets and the described methodology for inversion, that the operational monitoring of the aquifer can be done using A-DInSAR with ascending and/or descending satellite radar images. GNSS, using continuous or survey observation mode, can be used for validation and scaling purposes. GNSS stations should be installed in those locations that will best constrain the A-DInSAR results in areas of zero deformation for reference, maximum deformation areas for scaling and study of the time variation, or in low coherence areas so that both techniques complement each other. Continuous GNSS observations are preferable, if possible. Their proposed methodology would potentially reduce the cost of the geodetic monitoring system in a very important way. This effective A-DInSAR monitoring can be accomplished with the freely-available Copernicus Sentinel-1A and -1B satellite data, considering their global coverage and repeatability, ensuring their effective use for monitoring on a global scale [5].

4. General Summary

We have described a new inversion methodology, first presented by Camacho et al. [1], for modeling space and terrestrial geodetic displacement and gravity changes data in active areas with application tests. We assume that surface deformation and gravity changes are due to density and/or pressure changes allocated within extended source structures below the surface.

The characteristics of this approach can be summarized as follows:

(a): It permits simultaneous inversion of the 1D to 3D displacement data coming from different techniques, including terrestrial and space data (GNSS, DInSAR, leveling, etc.);
(b): Non-planar, non-gridded, inexact data can be used;
(c): It allows for objective modelling of two causative structures: Pressure and mass anomalies bodies;
(d): Simple analytical and well-known expressions are used for direct calculation for mass and pressure cells and their aggregation;
(e): This methodology works in a fully 3D context;
(f): There is not additional a priori requirements on the geometry and types of the causative sources;
(g): The method automatically determines the type and number, modeling the geometry, of the causative source structures;
(h): A free geometry of the causative structures is described by aggregation of small elemental cells. Surface deformation can be computed by adding the influence of the small considered prisms. Given that the used direct models are linear and the entire subsurface is assumed to be isotropic, superposition is allowable. Using this assumption, these linear equations permit the computation of surface deformation based on the superposition of many prismatic blocks within a compacting reservoir of any geometric shape [5]. The effect of each cell is computed using point sources [1,5].

When considering the validity of the results depending on the combination of the detected sources (if we have more than one pressure source), we must take into account the pressure change values employed and the adjusted distances between the different 3D pressure sources, which should be greater than 8 time their radii [124]. Solutions diverging from this criteria must be considering carefully and as a first, rough approach to the reality. Figure 26 show a schematic flow diagram of proposed inversion methodology.

5. Conclusions

In active areas, the source modelling for geodetic (terrestrial and space) data fitting is usually based on a priori fixed geometries (spherical, ellipsoidal prolate or oblate, etc.) for point or extended sources/causative bodies characterized by anomalous density or pressure. This new approach avoids this limitation to model the volume of these anomalous extended bodies in a 3D context with a free geometry.

In a first version of the methodology, the Mogi model was used as analytical solution for the nucleus of strain model. The Mogi model is a linear, first order solution. In a second version of the methodology, Geertsma’s nucleus of strain was implemented as the direct model solution, considering a poroelastic medium in place of a purely elastic one. This new version of the methodology allows us to consider the case of deformation associated with groundwater aquifer exploitation, increasing the applicability of the model/inversion methodology.

The described methodology is being implemented, PAF software [54], in the framework of the ESFRI (European Strategy Forum on Research Infrastructures) infrastructure EPOS [125] in the 3D-def service and will be available for a general use shortly.

Most of the active regions are characterized by complicated pattern of ground deformation resulting from multiple natural (e.g., volcanic loading producing inflation, deflation, dike intrusion, active faulting, flank instability and landslides) and anthropogenic (e.g., local uplift or subsidence) sources. We are developing an improvement of the described methodology which simultaneously estimates the typical ground deformation components in addition to the pressure source effects, allowing for a more general inverse methodology that includes dislocation sources, such as those given by the Okada’s expressions [126].

Author Contributions

J.F. and A.G.C. have contributed equally in designing and writing this manuscript.

Funding

This research was mainly funded by the Spanish Ministerio de Ciencia, Innovación y Universidades research project DEEP-MAPS, grant agreement number RTI2018-093874-B-I00, and also partially by the EU VII Framework Program, ESFRI, EPOS IP, grant agreement number: 676564-EPOS IP and the Spanish Ministerio de Economía y Competitividad Thematic Network EPOS Spain, grant number CGL2016-81965-REDT.

Acknowledgments

We thank the English revision carried out by K.F. Tiampo.

Conflicts of Interest

The authors declare no conflict of interest.

References

Camacho, A.G.; González, P.J.; Fernández, J.; Berrino, G. Simultaneous inversion of surface deformation and gravity changes by means of extended bodies with a free geometry: Application to deforming calderas. J. Geophys. Res. 2011, 116, B10401. [Google Scholar] [CrossRef]
Dzurisin, D. Volcano Deformation: Geodetic Monitoring Techniques; Springer-Praxis Books in Geophysical Sciences. Praxis Publishing Ltd.: Chichester, UK, 2007; p. 441. [Google Scholar]
Herring, T. Treatise on Geophysics, Geodesy, 1st ed.; Elsevier: Amsterdam, The Netherlands, 2009; Volume 3, p. 446. [Google Scholar]
Fernández, J.; Pepe, A.; Poland, M.P.; Sigmundsson, F. Volcano geodesy: Recent developments and future challenges. J. Volcanol. Geotherm. Res. 2017, 344, 1–12. [Google Scholar] [CrossRef]
Fernández, J.; Prieto, J.F.; Escayo, J.; Camacho, A.G.; Luzón, F.; Tiampo, K.F.; Palano, M.; Abajo, T.; Pérez, E.; Velasco, J.; et al. Modeling the two-and three-dimensional displacement field in Lorca, Spain, subsidence and the global implications. Sci. Rep. 2018, 8, 14782. [Google Scholar] [CrossRef] [PubMed]
Huang, Y.; Yu, M.; Xu, Q.; Sawada, K.; Moriguchi, S.; Yashima, A.; Liu, C.; Xue, L. InSAR-derived digital elevation models for terrain change analysis of earthquake-triggered flow-like kandslides bases on ALOS/PALSAR imagery. Environ. Earth Sci. 2015, 73, 7661–7668. [Google Scholar] [CrossRef]
Yu, M.; Huang, Y.; Zhou, J.; Mao, L. Modeling of landslide topography based on micro-unmanned aerial vehicle photography and structure-from-motion. Environ. Earth Sci. 2017, 76, 520. [Google Scholar] [CrossRef]
Biggs, J.; Ebmeier, S.K.; Aspinall, W.P.; Lu, Z.; Pritchard, M.E.; Sparks, R.S.J.; Mather, T.A. Global link between deformation and volcanic eruption quantified by satellite imagery. Nat. Commun. 2014, 5, 3471. [Google Scholar] [CrossRef] [PubMed]
Biggs, J.; Pritchard, M.E. Global volcano monitoring: What does it mean when volcanoes deform? Elements 2017, 13, 17–22. [Google Scholar] [CrossRef]
Tralli, D.M.; Blom, R.G.; Zlotnicki, V.; Donnellan, A.; Evans, D.L. Satellite remote sensing of earthquake, volcano, flood, landslide and coastal inundation hazards. ISPRS J. Photogramm. Remote Sens. 2005, 59, 185–198. [Google Scholar] [CrossRef]
Bilham, R. The seismic future of cities. Bull. Earthq. Eng. 2009, 7, 839–887. [Google Scholar] [CrossRef] [Green Version]
González, P.J. Medida y Caracterización de Deformaciones Usando Técnicas Geodésicas y de Teledetección. Aplicación en Volcanología y Sismotectónica. Ph.D. Thesis, Universidad Complutense de Madrid, Madrid, Spain, 2010. (In Spanish). [Google Scholar]
Fernández, J.; González, P.J.; Camacho, A.G.; Prieto, J.F.; Bru, G. An overview of geodetic volcano research in the Canary Islands. Pure Appl. Geophys. 2015, 172, 3189–3228. [Google Scholar] [CrossRef]
Terranova, C.; Ventura, G.; Vilardo, G. Multiple causes of ground deformation in the Napoli metropolitan area (Italy) from integrated persistent scatterers DinSAR, geological, hydrological, and urban infrastructure data. Earth-Sci. Rev. 2015, 146, 105–119. [Google Scholar] [CrossRef]
Sparks, R.S.J.; Biggs, J.; Neuberg, J.W. Monitoring volcanoes. Science 2012, 335, 1310–1311. [Google Scholar] [CrossRef] [PubMed]
Sigmundsson, F.; Hreinsdóttir, S.; Hooper, A.; Árnadóttir, T.; Pedersen, R.; Roberts, M.J.; Óskarsson, N.; Auriac, A.; Decriem, J.; Einarsson, P.; et al. Intrusion triggering of the 2010 Eyjafjallajökull explosive eruption. Nature 2010, 468, 426–430. [Google Scholar] [CrossRef] [PubMed]
Cannavò, F.; Camacho, A.G.; González, P.J.; Mattia, M.; Puglisi, G.; Fernández, J. Real time tracking of magmatic intrusions by means of ground deformation modeling during volcanic crises. Sci. Rep. 2015, 5, 10970. [Google Scholar] [CrossRef] [PubMed]
Galloway, D.; Jones, D.R.; Ingebritsen, S. Land Subsidence in the United States; US Geological Survey Circular: Reston, VA, USA, 1999. [Google Scholar]
Dokka, R.K. Modern-day tectonic subsidence in coastal Louisiana. Geology 2006, 34, 281–284. [Google Scholar] [CrossRef]
Amelung, F.; Galloway, D.L.; Bell, J.W.; Zebker, H.A.; Laczniak, R.J. Sensing the ups and downs of Las Vegas: InSAR reveals structural control of land subsidence and aquifer-system deformation. Geology 1999, 27, 483–486. [Google Scholar] [CrossRef]
González, P.J.; Fernández, J. Drought-driven transient aquifer compaction imaged using multitemporal satellite radar interferometry. Geology 2011, 39, 551–554. [Google Scholar] [CrossRef]
Boní, R.; Herrera, G.; Meisina, C.; Notti, D.; Béjar-Pizarro, M.; Zucca, F.; González, P.J.; Palano, M.; Tomás, R.; Fernández, J.; et al. Twenty-year advanced DInSAR analysis of severe land subsidence: The Alto Guadalentín basin (Spain) case study. Eng. Geol. 2015, 198, 40–52. [Google Scholar] [CrossRef]
Samsonov, S.V.; Tiampo, K.F.; Feng, W. Fast subsidence in downtown of Seattle observed with satellite radar. Remote Sens. Appl. Soc. Environ. 2016, 4, 179–187. [Google Scholar] [CrossRef]
Bru, G.; González, P.J.; Mateos, R.M.; Roldán, F.J.; Herrera, G.; Béjar-Pizarro, M.; Fernández, J. A-DInSAR monitoring of landslide and subsidence activity: A case of urban damage in arcos de la frontera, Spain. Remote Sens. 2017, 9, 787. [Google Scholar] [CrossRef]
Gourmelen, N.; Amelung, F.; Casu, F.; Manzo, M.; Lanari, R. Mining-related ground deformation in Crescent Valley, Nevada: Implications for sparse GPS networks. Geophys. Res. Lett. 2007, 34, L09309. [Google Scholar] [CrossRef]
Ma, C.; Cheng, X.; Yang, Y.; Zhang, X.; Guo, Z.; Zou, Y. Investigation on mining subsidence based on multi-temporal InSAR and time-series analysis of the small baseline subset—case study of working faces 22201-1/2 in Bu’ertai mine, shendong coalfield. Remote Sens. 2016, 8, 951. [Google Scholar] [CrossRef]
Lisowski, M. Analytical volcano deformation source parameters in volcano deformation: New geodetic monitoring techniques. In Volcano Deformation: Geodetic Monitoring Techniques; Dzurisin, D., Ed.; Springer-Verlag: Berlin/Heidelberg, Germany, 2007; pp. 279–304. [Google Scholar]
Segall, P. Earthquake and Volcano Deformation; Princeton University Press: Princeton, NJ, USA; Oxford, UK, 2010; p. 432. ISBN 978-0-691-13302-7. [Google Scholar]
Rymer, H.; Williams-Jones, G. Volcanic eruption prediction: Magma chamber physics from gravity and deformation measurements. Geophys. Res. Lett. 2000, 27, 16–2389. [Google Scholar] [CrossRef]
Fernández, J.; Tiampo, K.F.; Jentzsch, G.; Charco, M.; Rundle, J.B. Inflation or deflation? New results for Mayon volcano applying elastic-gravitational modeling. Geophys. Res. Lett. 2001, 28, 2349–2352. [Google Scholar] [CrossRef]
Masterlark, T. Magma intrusion and deformation predictions: Sensitivities to the Mogi assumptions. J. Geophys. Res. 2007, 112, B06419. [Google Scholar] [CrossRef]
Long, S.M.; Grosfils, E.B. Modeling the effect of layered volcanic material on magma reservoir failure and associated deformation, with application to Long Valley caldera, California. J. Volcanol. Geotherm. Res. 2009, 186, 349–360. [Google Scholar] [CrossRef]
Castaldo, R.; Gola, G.; Santilano, A.; De Novellis, V.; Pepe, S.; Manzo, M.; Manzella, A.; Tizzani, P. The role of thermo-rheological properties of the crust beneath Ischia island (Southern Italy) in the modulation of the ground deformation pattern. J. Volcanol. Geotherm. Res. 2017, 344, 154–173. [Google Scholar] [CrossRef]
Trasatti, E.; Giunchi, C.; Agostinetti, N.P. Numerical inversion of deformation caused by pressure sources: Application to Mount Etna (Italy). Geophys. J. Int. 2008, 172, 873–884. [Google Scholar] [CrossRef]
Battaglia, M.; Gottsmann, J.; Carbone, D.; Fernández, J. 4D volcano gravimetry. Geophysics 2008, 73, WA3–WA18. [Google Scholar] [CrossRef] [Green Version]
Carbone, D.; Poland, M.P.; Diament, M.; Grecco, F. The added value of 4D microgravimetry to the understanding of how volcanoes work. Earth Sci. Rev. 2017, 169, 146–179. [Google Scholar] [CrossRef]
Charco, M.; Fernández, J.; Luzón, F.; Rundle, J.B. On the relative importance of selfgravitation and elasticity in modeling volcanic ground deformation and gravity changes. J. Geophys. Res. 2006, 111, B03404. [Google Scholar] [CrossRef]
Johnson, D.J.; Eggers, A.A.; Bagnardi, M.; Battaglia, M.; Poland, M.P.; Miklius, A. Shallow magma accumulation at Kīlauea Volcano, Hawaii, revealed by microgravity surveys. Geology 2010, 38, 1139–1142. [Google Scholar] [CrossRef]
Bagnardi, M.; Poland, M.P.; Carbone, D.; Baker, S.; Battaglia, M.; Amelung, F. Gravity changes and deformation at Kīlauea Volcano, Hawaii, associated with summit eruptive activity, 2009–2012. J. Geophys. Res. 2014, 119, 7288–7305. [Google Scholar] [CrossRef]
Currenti, G. Numerical evidence enabling reconciliation of gravity and height changes in volcanic areas. Geophys. J. Int. 2014, 197, 164–173. [Google Scholar] [CrossRef]
Charco, M.; Camacho, A.G.; Tiampo, K.F.; Fernández, J. Spatiotemporal gravity changes on volcanoes: Assessing the importance of topography. Geophys. Res. Lett. 2009, 36, L08306. [Google Scholar] [CrossRef]
Fokker, P.A.; Wassing, B.B.T.; van Leijen, F.J.; Hanssen, R.F.; Nieuwland, D.A. Application of an ensemble smoother with multiple data assimilation to the Bergermeer gas field, using PS-InSAR. Geomech. Energy Environ. 2016, 5, 16–28. [Google Scholar] [CrossRef]
Vasco, D.W.; Wicks, C.; Karasaki, K.; Marques, O. Geodetic imaging: Reservoir monitoring using satellite interferometry. Geophys. J. Int. 2002, 149, 555–571. [Google Scholar] [CrossRef]
Tiampo, K.F.; Ouegnin, F.-A.; Valluri, S.; Samsonov, S.; Fernández, J.; Kapp, G. An elliptical model for deformation due to groundwater fluctuations. Pure Appl. Geophys. 2012, 169, 1443–1456. [Google Scholar] [CrossRef]
Geertsma, J. Land subsidence above compacting oil and gas reservoirs. J. Pet. Technol. 1973, 25, 734–744. [Google Scholar] [CrossRef]
Geertsma, J.; Van Opstal, G. A numerical technique for predicting subsidence above compacting reservoirs based on the nucleus of strain concept. Verh. Kon. Ned. Geol. Mijnbouwk 1973, 28, 63–78. [Google Scholar]
Segall, P. Stress and subsidence resulting from subsurface fluid withdrawal in the epicentral region of the 1983 Coalinga Earthquake. J. Geophys. Res. 1985, 90, 6801–6816. [Google Scholar] [CrossRef]
Vasco, D.W.; Karasaki, K.; Doughty, C. Using surface deformation to image reservoir dynamics. Geophysics 2000, 65, 132–147. [Google Scholar] [CrossRef]
Walsh, J.B. Subsidence above a planar reservoir. J. Geophys. Res. Solid Earth 2002, 107, 2202. [Google Scholar] [CrossRef]
Brown, N.J.; Woods, A.W.; Neufeld, J.A.; Richardson, C. Constraining Surface Deformation Predictions Resulting From Coal Seam Gas Extraction; Geoscience Australia: Symonston, Australia, 2014. [Google Scholar] [CrossRef]
González, P.J.; Tiampo, K.F.; Palano, M.; Cannavó, F.; Fernández, J. The 2011 Lorca earthquake slip distribution controlled by groundwater crustal unloading. Nat. Geosci. 2012, 5, 821–825. [Google Scholar] [CrossRef] [Green Version]
Charco, M.; Galán del Sastre, P. Efficient inversion of 3D finite element models of volcano deformation. Geophys. J. Int. 2014, 196, 1441–1454. [Google Scholar] [CrossRef]
Samsonov, S.V.; Tiampo, K.F.; Camacho, A.G.; Fernández, J.; González, P.J. Spatiotemporal analysis and interpretation of 1993–2013 ground deformation at Campi Flegrei, Italy, observed by advanced DInSAR. Geophys. Res. Lett. 2014, 41, 6101–6108. [Google Scholar] [CrossRef]
Camacho, A.G.; Fernández, J.; Cannavò, F. PAF: A software tool to estimate free-geometry extended bodies of anomalous pressure from surface deformation data. Comput. Geosci. 2018, 111, 235–243. [Google Scholar] [CrossRef]
Camacho, A.G.; Montesinos, F.G.; Vieira, R. A 3-D gravity inversion by means of growing bodies. Geophysics 2000, 65, 95–191. [Google Scholar] [CrossRef]
Davis, R.O.; Selvadurai, A.P.S. Elasticity and Geomechanics; Cambridge University Press: Cambridge, UK, 1996; p. 201. [Google Scholar]
Saleh, B. Underground deformation measurements using new quarts instruments. In Proceedings of the 95th Annual CIG Geomatics Conference, Can. Inst. of Geomatics, Ottawa, ON, Canada, 8–23 July 2002. [Google Scholar]
Rundle, J.B. Gravity changes and the Palmdale uplift. Geophys. Res. Lett. 1978, 5, 41–44. [Google Scholar] [CrossRef]
Walsh, J.B.; Rice, J.R. Local changes in gravity resulting from deformation. J. Geophys. Res. 1979, 84, 165–170. [Google Scholar] [CrossRef] [Green Version]
Currenti, G.; Del Negro, C.; Ganci, G. Modelling of ground deformation and gravity fields using finite element method: An application to Etna volcano. Geophys. J. Int. 2007, 169, 775–786. [Google Scholar] [CrossRef]
Williams, C.A.; Wadge, G. The effects of topography on magma chamber deformation models: Application to Mt. Etna and radar interferometry. Geophys. Res. Lett. 1998, 25, 1549–1552. [Google Scholar] [CrossRef]
Battaglia, M.; Hill, D.P. Analytical modeling of gravity changes and crustal deformation at volcanoes: The Long Valley caldera, California, case study. Tectonophysics 2009, 471, 45–57. [Google Scholar] [CrossRef]
Biot, M.A. General theory of three-dimensional consolidation. J. Appl. Phys. 1941, 12, 155–164. [Google Scholar] [CrossRef]
Tiede, C.; Tiampo, K.F.; Fernández, J.; Gerstenecker, C. Deeper understanding of non-linear geodetic data inversion using a quantitative sensitivity analysis. Nonlinear Processes Geophys. 2005, 12, 373–379. [Google Scholar] [CrossRef] [Green Version]
Farquharson, C.G.; Oldenbourg, D.W. Non-linear inversion using general measures of data misfit and model structure. Geophys. J. Int. 1998, 134, 213–227. [Google Scholar] [CrossRef] [Green Version]
Bertete-Aguirre, H.; Cherkaev, E.; Oristaglio, M. Non-smooth gravity problem with total variation penalization functional. Geophys. J. Int. 2002, 149, 499–507. [Google Scholar] [CrossRef] [Green Version]
Camacho, A.G.; Nunes, J.C.; Ortiz, E.; França, Z.; Vieira, R. Gravimetric determination of an intrusive complex under the Island of Faial (Azores): Some methodological improvements. Geophys. J. Int. 2007, 171, 478–494. [Google Scholar] [CrossRef]
Tarantola, A. Inverse Problem Theory; Elsevier: Amsterdam, The Netherlands, 1987; p. 613. [Google Scholar]
Camacho, A.G.; Montesinos, F.G.; Vieira, R. A 3-D gravity inversion tool based on exploration of model possibilities. Comput. Geosci. 2002, 28, 191–204. [Google Scholar] [CrossRef]
Gottsmann, J.; Camacho, A.G.; Marti, J.; Wooller, L.; Fernández, J.; Garcia, A.; Rymer, H. Shallow structure beneath the central volcanic complex of tenerife from new gravity data: Implications for its evolution and recent reactivation. Phys. Earth Planet. Inter. 2008, 168, 212–230. [Google Scholar] [CrossRef]
Corrado, G.; Guerra, I.; Lo Bascio, A.; Luongo, G.; Rampoldi, R. Inflation and microearthquake activity of Phlegraean Fields, Italy. Bull. Volcanol. 1976, 40, 1–20. [Google Scholar] [CrossRef]
Berrino, G.; Corrado, G.; Luongo, G.; Toro, B. Ground deformation and gravity changes accompanying the 1982 Pozzuoli Uplift. Bull. Volcanol. 1984, 47, 187–200. [Google Scholar] [CrossRef]
Orsi, G.; Civetta, L.; Del Gaudio, C.; De Vita, S.; Di Vito, M.A.; Isaia, R.; Petrazzuoli, S.M.; Ricciardi, G.P.; Ricco, C. Short-term ground deformation and seismicity in the resurgent Campi Flegrei caldera (Italy): An example of active block-resurgence in a densely populated area. J. Volcanol. Geotherm. Res. 1999, 91, 415–451. [Google Scholar] [CrossRef]
Gottsmann, J.; Berrino, G.; Rymer, H.; William-Jones, G. Hazard assessment during caldera unrest at the Campi Flegrei, Italy: A contribution from gravity-height gradients. Earth Planet. Sci. Lett. 2003, 211, 295–309. [Google Scholar] [CrossRef]
Berrino, G. Gravity changes induced by height-mass variations at the Campi Flegrei caldera. J. Volcanol. Geotherm. Res. 1994, 61, 293–309. [Google Scholar] [CrossRef]
Manconi, A.; Walter, T.R.; Manzo, M.; Zeni, G.; Tizzani, P.; Sansosti, E.; Lanari, R. On the effects of 3-D mechanical heterogeneities at Campi Flegrei caldera, southern Italy. J. Geophys. Res. 2010, 115, B08405. [Google Scholar] [CrossRef]
Lanari, R.; Berardino, P.; Borgström, S.; Del Gaudio, C.; De Martino, P.; Fornaro, G.; Guarino, S.; Ricciardi, G.P.; Sansosti, E.; Lundgren, P. The use of IFSAR and classical geodetic techniques for caldera unrest episodes: Application to the Campi Flegrei uplift event of 2000. J. Volcanol. Geotherm. Res. 2004, 133, 247–260. [Google Scholar] [CrossRef]
Casu, F.; Manzo, M.; Lanari, R. A quantitative assessment of the SBAS algorithm performance for surface deformation retrieval from DInSAR data. Remote Sens. Environ. 2006, 102, 195–210. [Google Scholar] [CrossRef]
Gottsmann, J.; Rymer, H.; Berrino, G. Unrest at the campi flegrei caldera (Italy): A critical evaluation of source parameters from geodetic data inversion. J. Volcanol. Geotherm. Res. 2006, 150, 132–145. [Google Scholar] [CrossRef]
Berrino, G.; Rymer, H.; Brown, G.C.; Corrado, G. Gravity height correlations for unrest calderas. J. Volcanol. Geotherm. Res. 1992, 53, 11–26. [Google Scholar] [CrossRef]
Talwani, P.; Acree, S. Pore pressure diffusion and the mechanism of reservoir-induced seismicity. Pure Appl. Geophys. 1984, 122, 947–965. [Google Scholar] [CrossRef]
Todesco, M.; Berrino, G. Modeling hydrothermal fluid circulation and gravity signals at the Phlegraean Fields caldera. Earth Planet. Sci. Lett. 2005, 240, 328–338. [Google Scholar] [CrossRef]
Berrino, G.; Corrado, G.; Riccardi, U. Sea gravity data in the Gulf of Naples: A contribution to delineating the Phlegraean Volcanic District. J. Volcanol. Geotherm. Res. 2008, 175, 241–252. [Google Scholar] [CrossRef]
Gottsmann, J.; Folch, A.; Rymer, H. Unrest at Campi Flegrei: A contribution to the magmatic versus hydrothermal debate from inverse and finite element modeling. J. Geophys. Res. 2006, 111, B07203. [Google Scholar] [CrossRef]
Battaglia, M.; Obrizzo, F.; Pingue, F.; De Natale, G. Evidence for fluid migration as the source of deformation at Campi Flegrei caldera (Italy). Geophys. Res. Lett. 2006, 33, L01307. [Google Scholar] [CrossRef]
Berardino, P.; Fornaro, G.; Lanari, R. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Trans. Geosci. Rem. Sens. 2002, 40, 2375–2383. [Google Scholar] [CrossRef] [Green Version]
Usai, S. A least squares database approach for SAR interferometric data. IEEE Trans. Geosci. Rem. Sens. 2003, 41, 753–760. [Google Scholar] [CrossRef] [Green Version]
Samsonov, S.; van der Koij, M.; Tiampo, K. A simultaneous inversion for deformation rates and topographic errors of DInSAR data utilizing linear least square inversion technique. Comput. Geosci. 2011, 37, 1083–1091. [Google Scholar] [CrossRef]
Ferretti, A.; Prati, C.; Rocca, F. Permanent scatterers in SAR interferometry. IEEE Trans. Geosci. Rem. Sens. 2001, 39, 8–20. [Google Scholar] [CrossRef]
Hooper, A. A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches. Geophys. Res. Lett. 2008, 35, L16302. [Google Scholar] [CrossRef]
Samsonov, S.; d’Oreye, N. Multidimensional time series analysis of ground deformation from multiple InSAR data sets applied to Virunga volcanic province. Geophys. J. Int. 2012, 191, 1095–1108. [Google Scholar] [CrossRef]
Samsonov, S.; d’Oreye, N.; Smets, B. Ground deformation associated with post-mining activity at the French-German border revealed by novel InSAR time series method. Int. J. Appl. Earth Obs. Geoinf. 2013, 23, 142–154. [Google Scholar] [CrossRef]
Samsonov, S.; d’Oreye, N.; González, P.; Tiampo, K.; Ertolahti, L.; Clague, J.J. Rapidly accelerating subsidence in the Greater Vancouver region from two decades of ERS-ENVISAT-RADARSAT-2 DInSAR measurements. Rem. Sens. Environ. 2014, 143, 180–191. [Google Scholar] [CrossRef] [Green Version]
Amoruso, A.; Crescentini, L.; Sabbetta, I. Paired deformation sources of the Campi Flegrei caldera (Italy) required by recent (1980–2010) deformation history. J. Geophys. Res. Solid Earth 2014, 119, 858–879. [Google Scholar] [CrossRef]
Chiodini, G.; Caliro, S.; De Martino, P.; Avino, R.; Gherardi, F. Early signals of new volcanic unrest at Campi Flegrei caldera? Insights from geochemical data and physical simulations. Geology 2012, 40, 943–946. [Google Scholar] [CrossRef]
Aloisi, M.; Bonaccorso, A.; Cannavò, F.; Gambino, S.; Mattia, M.; Puglisi, G.; Boschi, E. A new dyke intrusion style for the Mount Etna May 2008 eruption modelled through continuous tilt and GPS data. Terr. Nova 2009, 21, 316–321. [Google Scholar] [CrossRef]
Palano, M.; Rossi, M.; Cannavò, F.; Bruno, V.; Marco, A.; Daniele, P.; Mario, P.; Siligato, G.; Mattia, M. Etn@ref: A geodetic reference frame for Mt. Etna GPS networks. Ann. Geophys. 2010, 53, 49–57. [Google Scholar]
Bonaccorso, A.; Bonforte, A.; Currenti, G.; Del Negro, C.; Di Stefano, A.; Greco, F. Magma storage, eruptive activity and flank instability: Inferences from ground deformation and gravity changes during the 1993-2000 recharging of Mt. Etna volcano. J. Volcanol. Geotherm. Res. 2011, 200, 245–254. [Google Scholar] [CrossRef]
Maronna, R.A.; Martin, D.R.; Yohai, V.J. Robust Statistics: Theory and Methods; John Wiley & Sons Ltd.: Chichester, UK, 2006; p. 436. [Google Scholar]
Patané, D.; De Gori, P.; Chiarabba, C.; Bonaccorso, A. Magma ascent and the pressurization of Mount Etna’s volcanic system. Science 2003, 290, 2061–2063. [Google Scholar] [CrossRef]
Aloisi, M.; Mattia, M.; Monaco, C.; Pulvirenti, F. Magma, faults, and gravitational loading at Mount Etna: The 2002–2003 eruptive period. J. Geophys. Res. 2011, 116, B05203. [Google Scholar] [CrossRef]
D’Auria, L.; Giudicepietro, F.; Martini, M.; Lanari, R. The 4D imaging of the source of ground deformation at Campi Flegrei caldera (southern Italy). J. Geophys. Res. 2012, 117, B08209. [Google Scholar] [CrossRef]
Fernández, J.; Rundle, J.B. Gravity changes and deformation due to a magmatic intrusion in a two-layered crustal model. J. Geophys. Res. 1994, 99, 2737–2746. [Google Scholar] [CrossRef] [Green Version]
Rundle, J.B. Viscoelastic-gravitational deformation by a rectangular thrust fault in a layered Earth. J. Geophys. Res. 1982, 87, 7787–7796. [Google Scholar] [CrossRef]
AEMET, Agencia Estatal de Meteorología. Available online: www.aemet.es/es/datos_abiertos/ AEMET_OpenData (accessed on 15 January 2018).
Bourgois, J.; Mauffret, A.; Ammar, A.; Demnati, A. Multichannel seismic data imaging of inversion tectonics of the Alboran Ridge (western Mediterranean Sea). Geo Marine Lett. 1992, 12, 117–122. [Google Scholar] [CrossRef]
Martínez-Díaz, J.J. Stress field variation related to fault interaction in a reverse oblique-slip fault: The Alhama de Murcia fault, Betic Cordillera, Spain. Tectonophysics 2002, 356, 291–305. [Google Scholar] [CrossRef]
Palano, M.; González, P.J.; Fernández, J. Strain and stress fields along the Gibraltar Orogenic Arc: Constraints on active geodynamics. Gondwana Res. 2013, 23, 1071–1088. [Google Scholar] [CrossRef]
Cerón, J.C.; Pulido-Bosch, A. Groundwater problems resulting from CO 2 pollution and overexploitation in Alto Guadalentín aquifer (Murcia, Spain). Environ. Geol. 1996, 28, 223–228. [Google Scholar] [CrossRef]
Confederación Hidrográfica del Segura. Plan especial de actuación en situaciones de alerta y eventual Sequía. Tech. Rep. 2006, 298. [Google Scholar]
IGN. Instituto Geográfico Nacional. Available online: http://www.ign.es/ (accessed on 15 January 2018).
Tomás, R.; Herrera, G.; Lopez-Sanchez, J.; Vicente-Guijalba, F.; Cuenca, A.; Mallorqui, J.J. Study of the land subsidence in Orihuela City (SE Spain) using PSI data: Distribution, evolution and correlation with conditioning and triggering factors. Eng. Geol. 2010, 115, 105–121. [Google Scholar] [CrossRef]
Herrera, G.; Tomás, R.; Monells, D.; Centolanza, G.; Mallorquí, J.J.; Vicente, F.; Navarro, V.D.; Lopez-Sanchez, J.M.; Sanabria, M.; Cano, M.; et al. Analysis of subsidence using TerraSAR-X data: Murcia case study. Eng. Geol. 2010, 116, 284–295. [Google Scholar] [CrossRef]
Zhu, L.; Gong, H.; Li, X.; Wang, R.; Chen, B.; Dai, Z.; Teatini, P. Land subsidence due to groundwater withdrawal in the northern Beijing plain, China. Eng. Geol. 2015, 193, 243–255. [Google Scholar] [CrossRef]
Pacheco-Martínez, J.; Cabral-Cano, E.; Wdowinski, S.; Hernández-Marín, M.; Ortiz-Lozano, J.A.; Zermeño-de-León, M.E. Application of InSAR and gravimetry for land subsidence hazard zoning in Aguascalientes, Mexico. Remote Sens. 2015, 7, 17035–17050. [Google Scholar] [CrossRef]
Solari, L.; Ciampalini, A.; Raspini, F.; Bianchini, S.; Moretti, S. PSInSAR analysis in the Pisa urban area (Italy): A case study of subsidence related to stratigraphical factors and urbanization. Remote Sens. 2016, 8, 120. [Google Scholar] [CrossRef]
Chaussard, E.; Milillo, P.; Bürgmann, R.; Perissin, D.; Fielding, E.J.; Baker, B. Remote sensing of ground deformation for monitoring groundwater management practices: Application to the Santa Clara Valley during the 2012–2015 California drought. J. Geophys. Res. Solid Earth 2017, 122, 8566–8582. [Google Scholar] [CrossRef]
Chen, G.; Zhang, Y.; Zeng, R.; Yang, Z.; Chen, X.; Zhao, F.; Meng, X. Detection of land subsidence associated with land creation and rapid urbanization in the Chinese Loess Plateau using time series InSAR: A case study of Lanzhou New district. Remote Sens. 2018, 10, 270. [Google Scholar] [CrossRef]
Hanssen, R.F. Radar Interferometry: Data Interpretation and Error Analysis; Kluwer Academic Publishers: Dordrecht, The Netherlands, 2001; p. 308. [Google Scholar] [CrossRef]
Samieie-Esfahany, S.; Hanssen, R.F.; van Thienen-Visser, K.; Muntendam-Bos, A.; Systems, S. On the effect of horizontal deformation on InSAR subsidence estimates. In Proceedings of the Fringe Workshop, Frascati, Italy, 30 November–4 December 2009. [Google Scholar]
Mora, O.; Ordoqui, P.; Iglesias, R.; Blanco, P. Earthquake rapid mapping using ascending and descending Sentinel-1 TOPSAR interferograms. Procedia Comput. Sci. 2016, 100, 1135–1140. [Google Scholar] [CrossRef]
Blanco-Sánchez, P.; Mallorquí, J.J.; Duque, S.; Monells, D. The coherent pixels technique (CPT): An advanced DInSAR technique for nonlinear deformation monitoring. Pure Appl. Geophys. 2008, 165, 1167–1193. [Google Scholar] [CrossRef]
Ezquerro, P.; Guardiola-Albert, C.; Herrera, G.; Fernández-Merodo, J.A.; Béjar-Pizarro, M.; Bonì, R. Groundwater and subsidence modeling combining geological and multi-satellite SAR data over the alto guadalentín aquifer (SE Spain). Geofluids 2017, 1–17. [Google Scholar] [CrossRef]
Pascal, K.; Neuberg, J.; Rivalta, I. On precisely modelling surface deformation due to interacting magma chambers and dykes. Geophys. J. Int. 2014, 196, 253–278. [Google Scholar] [CrossRef]
European Plate Observing System. Available online: https://www.epos-ip.org/ (accessed on 25 July 2019).
Okada, Y. Surface deformation due to shear and tensile faults in a halfspace. Bull. Seismol. Soc. Amer. 1985, 75, 1135–1154. [Google Scholar]

Figure 1. (a) Location of survey benchmarks for repeated gravity and leveling in Campi Flegrei. (b) Temporal changes for gravity (red) and elevation (blue) at Serapeo from 1980 to 2000. The vertical dimension of the symbols is representative of the errors. A high correlation is observed between both data types, but the elevation values show a more continuous pattern. Modified from [1].

Figure 2. (a) LOS deformation velocity computed from ascending passes for the period 1993–2000 and (b) LOS deformation velocity computed from descending passes for the period 1992–2000. Modified from [1].

Figure 3. Three cross sections of the 3-D model for depressurization: Horizontal (depth 1500 m) and NNE–SSW and WNW–ESE vertical sections of the under-pressure model across a central position resulting from the simultaneous inversion of the gravity changes, leveling changes, and DInSAR data in Campi Flegrei for 1992–2000 assuming an elastic half-space [1].

Figure 4. Some additional views of the LOS ascending data fit: (a) residual map and (b) WE central profile for observed-modeled comparison [1].

Figure 5. MSBAS results, 1993–2013, for the images detailed in Table 1. (a) Vertical cumulative component of deformation in centimeters, 1993–2013; (b) east-west cumulative component of deformation in centimeters, 1993–2013; (c) time series of vertical and east-west components shown in panels (a) and (b) at location of maximum subsidence, identified with the green dot. The reference location for MSBAS processing is located at 4532380, 4335351 (14.23°N, 40.94°E). Modified from [53].

Figure 6. Source location, depth and shape for deflation period of 1993–1999. (a) 3D perspective; (b) map view of source below Campi Flegrei caldera; (c) EW vertical profile; (d) NS vertical profile [53].

Figure 7. (a) Observed vertical displacement rate (left), cm/yr, for the subsidence period, 1993–1999; modelled vertical displacement rate from inversion (center); and residual of observed and modelled displacements (right). (b) Observed EW displacement rate (left), cm/yr, 1993–1999; modelled EW displacement from inversion (center); and residual of observed and modelled displacements (right). Images were saturated at the corresponding scales in order to improve comparison and highlight differences among panels [53].

Figure 8. Source location, depth and shape for deflation period of 1999–2000. (a) 3D perspective; (b) map view of source below Campi Flegrei caldera; (c) EW vertical profile; (d) NS vertical profile [53].

Figure 9. (a) Observed vertical displacement rate (left), cm/yr, for the subsidence period, 1999–2000; modelled vertical displacement rate from inversion (center); and residual of observed and modelled displacements (right). (b) Observed EW displacement rate (left), cm/yr, 1999–2000; modelled EW displacement from inversion (center); and residual of observed and modelled displacements (right). Images were saturated at the corresponding scales in order to improve comparison and highlight differences among panels [53].

Figure 10. Source location, depth and shape for deflation period of 2000–2005. (a) 3D perspective; (b) map view of source below Campi Flegrei caldera; (c) EW vertical profile; (d) NS vertical profile [53].

Figure 11. (a) Observed vertical displacement rate (left), cm/yr, for the subsidence period, 2000–2005; modelled vertical displacement rate from inversion (center); and residual of observed and modelled displacements (right). (b) Observed EW displacement rate (left), cm/yr, 2000–2005; modelled EW displacement from inversion (center); and residual of observed and modelled displacements (right). Images were saturated at the corresponding scales in order to improve comparison and highlight differences among panels [53].

Figure 12. Source location, depth and shape for deflation period of 2005–2007. (a) 3D perspective; (b) map view of source below Campi Flegrei caldera; (c) EW vertical profile; (d) NS vertical profile [53].

Figure 13. (a) Observed vertical displacement rate (left), cm/yr, for the subsidence period, 2005–2007; modelled vertical displacement rate from inversion (center); and residual of observed and modelled displacements (right). (b) Observed EW displacement rate (left), cm/yr, 2005–2007; modelled EW displacement from inversion (center); and residual of observed and modelled displacements (right). Images were saturated at the corresponding scales in order to improve comparison and highlight differences among panels [53].

Figure 14. Source location, depth and shape for deflation period of 2007–2013. (a) 3D perspective; (b) map view of source below Campi Flegrei caldera; (c) EW vertical profile; (d) NS vertical profile [53].

Figure 15. (a) Observed vertical displacement rate (left), cm/yr, for the subsidence period, 2007–2013; modelled vertical displacement rate from inversion (center); and residual of observed and modelled displacements (right). (b) Observed EW displacement rate (left), cm/yr, 2007–2013; modelled EW displacement from inversion (center); and residual of observed and modelled displacements (right). Images were saturated at the corresponding scales in order to improve comparison and highlight differences among panels [53].

Figure 16. Planar view (a) and EW vertical cut view (b) of the main elements from the real time modelling process covering the inflation period and the eruption on 13 May 2008, and possible connections between the source bodies and the location of the earthquakes (blue circles) for 12–14 May 2008. Contours correspond to surface topography. Gray triangles indicate the location of the GPS stations. Shaded gray area outlines the position of the High Velocity Body (HVB) [100]. In green, the cumulated sources during the yearlong recharging phase and in orange the sources active during the day of the eruption. Cannavò et al. [17] suggested fast magma ascent from the reservoir level (2–3 km depth bsl) along paths indicated by the source geometries and the earthquake locations (dashed lines in the figures) [17].

Figure 17. (a) Best windowing. In blue the mean (for all the stations and components) standard deviations, for different window sizes, of the estimated displacement time series from GPS solutions. (b) Mean autocorrelation function with shaded error bar for all the 1-Hz GPS time series. It is possible to see the steep decay of autocorrelation after lags of a few hundreds of seconds. Modified from [17].

Figure 18. Model for the inflation phase. Cumulated magmatic body (in green) modelled for the period 1 January 2007 to 12 May 2008 and time evolution (in red) for the sub-period that corresponds to the months of June and July 2007, when the system was subject to a recharging episode, which appears as an ascending high-pressure body. Parallelepiped sources represent active point pressure sources within the same volume. Blue circles correspond to location of the earthquakes for 12–14 May 2008 [17].

Figure 19. Time sequence of pressurized source models (described by aggregation of parallelepiped cells which map active point pressure sources) corresponding to key instants from the real time inversions during the 13 May 2008, eruption of Mount Etna volcano at different UTC times: (a) at 7:00, (b) at 8:30, (c) at 11:30, and (d) at 15:00. Red volumes are positive-pressure sources while blue volumes are negative pressure sources. Structures marked in orange correspond to the cumulated sources describing the main source structures across the sequence. Contours correspond to surface topography. Gray triangles indicate the location of the GPS stations. Blue circles correspond to location of the earthquakes for 12–14 May 2008 [17].

Figure 20. Graphic summary of the sources modeled in the studied 2008 Mt. Etna eruption. HVB: High velocity body, F: the dyke as modelled by the present study and previous works [102], M: main source body for the inflation phase preceding the eruption. UTM coordinates and depth are expressed in meters. The azimuth view is 40° [17].

Figure 21. 3D uncertainty maps for the considered GPS network of 10 stations, representing, respectively, (a) the pressure changes, (b), the depth changes, and (c) the horizontal deviations required to produce a surface deformation with quadratic mean value 1 mm. Map (d) represents the minimum horizontal deviation for an isolated body of 1 km³ with a pressure change given in panel (a) to produce a deformation at the GPS stations with rms magnitude of 1 mm [17].

Figure 22. (a) Location of the Alto Guadalentín Basin, the Bajo Guadalentín Basin and the Guadalentín River that formed the two basins. Black lines depict main faults in the area. The locations and names of the main cities in the area are shown. The topography has been obtained from MDT05 2015 CC-BY 4.0 digital elevation model [111]. (b) Subsidence area detected in previous studies31 by means of DInSAR techniques along the Alto Guadalentín Basin. Subsidence rates have a maximum of 16 cm/yr for the period 2006–2011 located ~4 km south-west the city of Lorca. The black stars are damage locations due to the M = 5.1 May 2008 Lorca earthquake. Red lines are main faults (AMF, Alhama de Murcia Fault). The contour lines indicate 2 cm/yr DInSAR subsidence due to groundwater pumping. (c) Location of the monitoring GNSS control stations deployed in the area of Alto Guadalentín. The network consists of 33 monitoring stations (blue circles show their location) and covers an area of about 70 km². The network is designed to allow high accuracy GNSS surveys and also includes two existing continuous GNSS stations. Main population centers are depicted with white stars. Modified from [5].

Figure 23. Displacement rates determined from GNSS observations. Results corresponding to the period November 2015–February 2017. (a) Annual vertical displacement rates, subsidence, measured with standard confidence bars. (b) Average annual horizontal displacements with standard confidence regions. Additional results are shown in the Supplementary Information from [5]. (c) and (d) show the results obtained from the A-DInSAR processing using CPT technique. Both geometries, ascending and descending, have been processed using a multilook window of 3 × 13 pixels (azimuth × range) which generates a square pixel of about 60 × 60 meters in ground resolution. Coherence method has been used for pixel selection coherence method. Results are shown for the period November 2015–February 2017. (c) Line of Sight (LOS) velocity values obtained for the ascending orbit. (d) LOS velocity values for the descending orbit. Black dots locate the GNSS stations. Modified from [5].

Figure 24. East-West and Vertical displacements obtained by A-DInSAR. (a) Horizontal (East-West) and (b) vertical (Up-Down) displacement rates estimations obtained by decomposition of the LOS detected velocity using ascending and descending orbits. GNSS displacements are also plotted with arrows to compare. Results are shown for the period November 2015–February 2017. Modified from [5].

Figure 25. Representation of the inversion results obtained for the 1D, 2D, 3D and 2D+3D considered data sets. (a) Obtained source for Case A; (b) for Case B; (c) for Case C; (d) Results for Case D; (e) for Case E; (f) for Case F; (g) for Case G; (i) for Case I and (j) for Case J. Blue color indicates negative pressure value cells, produced by water extraction. White color indicates positive pressure change cells. These positive pressure sources adjust the errors and the effects of other deformation sources, different from water extraction (e.g., of tectonic origin) [5].

Figure 26. Schematic flow diagram of the described inversion methodology [1,5,17,53,54] where we start from the data set, the medium characteristics and its 3D gridding, to get (using the direct model equations and complementary conditions) a 3D source model of the anomalous sources via a growth process.

Table 1. Seven DInSAR data sets used in this study with continuous coverage from 1993 until 2013 ^a [53].

DInSAR Data Set	Orbit	Coverage Period ^b	θ	Φ	N	M
ERS, track 129	asc	10/1/1993–17/9/2008	344.1	23.2	90	215
ERS, track 036	dsc	8/6/1992–25/12/2008	194.1	23.2	84	164
Envisat, track 129	asc	13/11/2002–16/12/2009	344.0	22.8	46	120
Envisat, track 036	dsc	5/6/2003–21/10/2010	195.9	22.8	60	158
R2, S3	asc	19/1/2009–2/8/2013	348.7	35.1	42	166
R2, S3	dsc	27/12/2008–3/8/2013	190.4	35.1	53	290
R2, F6	asc	29/12/2008–5/8/2013	351.0	48.3	50	158

^a Incidence angle φ (degrees), azimuth angle θ (degrees), number of available single-look complex SAR images in each set N, and number of computed highly coherent interferograms M. Combined coverage: 10 January 1993 to 3 August 2013. Total number of unique time steps (measured daily) = 385 (48 repeated by different sensors). Total M = 1271. Ascending (asc) and descending (dsc). ^b Dates are formatted as day/month/year.

Table 2. Six time periods used for modeling expansion/contraction sources below Campi Flegrei Caldera and associated location, pressure, volume, and depth to the center of mass ^a [53].

Time Period (Years)	X Location UTM (m)	Y Location UTM (m)	Depth (m bsl)	Pressure (MPa × 10⁹/yr)	Volume (km³/yr)
	(±50)	(±50)	(±60)	(±0.05)	(±0.7)
1993–1999	426573	4519313	1372	−0.37	6.2
1999–2000	426419	4519781	1434	+0.51	8.5
2000–2005	426469	4519501	1874	−0.13	2.2
2005–2007	427388	4519326	2126	+0.69	11.4
2007–2007	426184	4519238	1683	−0.47	7.9
2007–2013	426669	4519315	1819	+0.49	8.2

^a UTM = Universal Transverse Mercator.

Table 3. Numerical summary of the inversion results for the ten cases considered ¹ [5].

CASE	Intensity MPa×Km³	Misfit (cm)	Mean Model Intensity(MPa × Km³)	Pres. (MPa)	Vol. (Km³)	Displacement Components Considered	Nº of Data Used
A	−41	0.36	−42.7±3.8 (9%) [−40±1.4 (4%)]	−3	14.2 [13.3]	1D	1505
B	−39	0.31					1203
C	−48	0.19					1572
D	−32	0.30	−33.9 ± 1.9 (6%) [−32.5 ± 1.1 (3%)]	−3	11.3 [10.8]	2D	1505
E	−31	0.28					1203
F	−33	0.32					2708
G	−37	0.45					3144
H	−34	0.94				3D	108
I	−34	0.43				2D + 3D	2816
J	−36	0.63				2D + 3D	3252

¹ Data are grouped into two types: 1D data (Cases A–C), 2D and 2D +3D (Cases D–J). 2D data come from A-DInSAR results. 2D+3D data sets combine data obtained from the A-DInSAR study with those obtained from GNSS observation campaigns. For each of these two sets, which give rise to very similar results, mean values of intensity and values of volume variation in the aquifer, as a function of the pressure variation assumed, are given. Figures in square brackets denote average values determined not considering Cases C, and Cases G, H and J respectively. See text for more details.

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Camacho, A.G.; Fernández, J. Modeling 3D Free-geometry Volumetric Sources Associated to Geological and Anthropogenic Hazards from Space and Terrestrial Geodetic Data. Remote Sens. 2019, 11, 2042. https://doi.org/10.3390/rs11172042

AMA Style

Camacho AG, Fernández J. Modeling 3D Free-geometry Volumetric Sources Associated to Geological and Anthropogenic Hazards from Space and Terrestrial Geodetic Data. Remote Sensing. 2019; 11(17):2042. https://doi.org/10.3390/rs11172042

Chicago/Turabian Style

Camacho, Antonio G., and José Fernández. 2019. "Modeling 3D Free-geometry Volumetric Sources Associated to Geological and Anthropogenic Hazards from Space and Terrestrial Geodetic Data" Remote Sensing 11, no. 17: 2042. https://doi.org/10.3390/rs11172042

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Modeling 3D Free-geometry Volumetric Sources Associated to Geological and Anthropogenic Hazards from Space and Terrestrial Geodetic Data

Abstract

1. Introduction

2. Inversion Approach for Deformation and Gravity Changes

2.1. System of Nonlinear Equations

2.2. Misfit Conditions

2.3. Exploration Approach for Solving the System

3. Application Results

3.1. Modeling of Campi Flegrei Unrest 1992–2000 Using Deformation and Gravity Changes

Discussion

3.2. Modeling of Campi Flegrei Unrest 1993–2013 Using Only Displacement Data

Discussion

3.3. Real Time Tracking of Magmatic Intrusions During Volcanic Crisis: 2008 May Etna Eruption

3.3.1. Results

3.3.2. Discussion

3.4. Land Subsidence Associated with Overexploitation of Aquifers: Lorca, Spain

Discussion and Inversion of the 3D Displacement Field

4. General Summary

5. Conclusions

Author Contributions

Funding

Acknowledgments

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI