Introduction

The modeling of the spatial and temporal evolution of the groud displacement field is relevant for the understanding of natural phenomena characterizing the Earth’s surface1,2,3,4,5,6. Parametric optimization, tomographic inversion and imaging methods are different approaches aimed at evaluating the physical and geometrical parameters of displacement sources, such as the ruptures of fault planes in case of earthquakes or the pressure/volume variations inside plumbing systems for volcanic unrests7,8,9,10,11,12,13. However, modeling strategies can fail in determining a unique solution of the source parameters when model constraints are not available, as can occur for volcanic unrest. Spatially dense ground displacement data along the three directions (3-D), i.e., East-West (E-W), North-South (N-S) and vertical (Up-Down, U-D), can undoubtedly help improve the modeling performance and evaluate results with more significant reliability14,15,16,17,18,19.

In this context, Synthetic Aperture Radar Interferometry (InSAR) has been playing a pivotal role in studying the Earth’s surface dynamics, allowing the detection and monitoring of surface movements20,21,22. The InSAR technology has extensively been applied with data collected by single SAR instruments, allowing the detection of the projection of the ground displacements along the radar-to-target line of sight (LOS) direction. Moreover, when at least two independent SAR datasets collected from multi-viewing, complementary illumination angles (i.e., from ascending and descending orbits) are available, the LOS-projected (1-D) ground displacements can be combined to recover estimates of the E-W and U-D components (2-D), with uncertainties of a few millimeters. The recent advances in SAR technology and the availability of extensive archives of SAR data collected by several satellite radar instruments have promoted the development of new approaches for the retrieval of 3-D ground displacements, using the complementary (amplitude and phase) pieces of information embedded in the different SAR datasets [e.g., 23,24,25]. However, because almost all modern spaceborne SAR satellites fly along near-polar orbits, the LOS-projected ground displacements are less sensitive to N-S displacements. To overcome this limitation, new technological solutions have been proposed, such as Pixel-Offset Tracking (POT)15,26, the Multiple Aperture Interferometry (MAI) techniques27,28 and the Burst Overlapped Interferometry (BOI)4,29,30. It has been shown that POT measurements have uncertainties of about 1/30th of the pixel size (e.g., ~ 10 cm for the COSMO-SkyMed and TerraSAR-X SAR systems). Conversely, MAI azimuthal displacement measurement precisions depend on the azimuthal bandwidth, the Effective Number of Looks (ENL), and the coherence of the imaged targets. For instance, for Sentinel-1 SAR data, the azimuthal measurement precision is about 3–4 cm with ENL = 20 and coherence values of 0.614,31. More recently, the enhanced spectral diversity and the BOI techniques have also been developed for the fine co-registration of SAR data pairs and the retrieval of the azimuthal ground displacements. In this case, BOI azimuthal ground displacements, albeit restricted only to the overlapped regions between consecutive bursts, have enhanced precisions of about 1.5 cm (for coherence values of 0.6 and ENL = 40). Nevertheless, for volcanic applications, the N-S ground displacement precisions and/or the covered region guaranteed by these technologies are hardly sufficient for a comprehensive understanding of phenomena occurring during the early pre-eruptive phases and, for instance, in these cases, the application of MAI, POT and BOI seems somewhat limited29.

In this work, we propose a new method for providing the 3-D ground displacement field that relies on LOS-projected InSAR measurements and applying potential field theory. Specifically, multi-viewing InSAR measurements are combined to recover precise (up to a few millimeters) ground displacement of the E-W and U-D components over a set of coherent targets. Then, the potential theory is used to infer the N-S component from the InSAR-driven E-W and U-D ones. This approach indeed leads to reliable imaging of the N-S displacement when the ground displacement is a potential field, namely for the approximated conditions of the elasticity theory32,33,34,35. The employment of the potential theory represents the innovative feature of the proposed approach since the imaging of the displacements along the N-S direction is performed only based on physical requirements characterizing the deformative mechanism; it is therefore not necessary to recur to independent measurements of the N-S components, for instance, using advanced SAR technologies.

We focus on volcano deformations and treat theoretical aspects to explain how to employ the potential theory for imaging the ground displacement components related to volcanic unrests. Application of potential theory to simulated noisy ground displacement datasets is here used to confirm the soundness of the described argumentations; specifically, we start from the analysis of the favorable scenario of hydrostatically pressurized sources with any geometry embedded in layered media to then address not-hydrostatically pressurized and multi-source cases for which the potential theory turns out to work only approximately. Preliminary analyses quantifying errors and uncertainties are also performed.

Finally, we used the designed procedure to infer the 3-D ground displacement field of Sierra Negra volcano. We considered the unrest event that started in February 2017, leading to the 26th June 2018 eruption, and we computed the N-S ground displacement velocity maps of the investigated area and the related time-series. We compared the estimated E-W and U-D and the retrieved N-S component’s displacement time series to available Global Navigation Satellite System (GNSS) data.

Methods

This section first summarizes the basic rationale of the InSAR technology adopted for retrieving the E-W and U-D ground displacement time-series. Then, we describe the potential field theory for elastic displacements to retrieve the N-S component. We also show applications of potential theory to simulated displacement fields in order to highlight its performance and uncertainties. Finally, both techniques are employed to reconstruct the 3-D ground displacements.

InSAR data processing

Ground displacements of the imaged areas are retrieved by applying the Differential Interferometric Synthetic Aperture Radar (DInSAR) technique21,36,37 to independent sets of Sentinel-1 SAR images, collected through complementary illumination angles (i.e., through ascending and descending flight tracks). DInSAR methods rely on the computation of the phase difference between pairs of temporally separated SAR images and allow retrieving a 1-D component of the ground displacement, representing the displacement projection along the radar-to-target LOS direction. For every single SAR dataset, the LOS-projected ground displacement time-series are recovered by applying the multi-temporal Small Baseline Subset (SBAS) technique38 to a set of short baselines, multi-looked (average) interferograms. For the interferogram generation, we employ precise orbital information and a digital elevation model (DEM) of the study area, e.g., the NASA Shuttle Radar Topography Mission (SRTM) DEM. Also, fine co-registration procedures39,40,41 are performed to process Sentinel-1 SAR data collected through the Terrain Observation with Progressive Scans (TOPS) mode. The generated interferograms are unwrapped42 and then inverted in a Least-Square (LS) sense through the SBAS method. Finally, the relevant time-series of the E-W and U-D ground displacement components are estimated for every coherent target seen from both (ascending/descending) flight paths by combining the observed LOS-projected measurements25,36,43 and profiting from the inherently different multi-platform acquisition geometries.

Potential theory for elastic displacement field

Let us consider a cartesian reference system, oriented so that the x-, y- and z-axes correspond to the East, North and vertical directions. The displacement vector field d(x,y,z) defines how any object is deformed, measuring the change of position of any of its points from an initial (xi,yi,zi) to a final (xf, yf, zf) set of coordinates. The x-, y- and z- differences between the final and initial positions describe the d(x,y,z) components, namely u(x,y,z) and v(x,y,z) for position change (i.e., displacement) along the East (i.e., E-W component) and North (i.e., N-S component) directions, and w(x,y,z) for altitude/depth changes (i.e., U-D component).

In order to show the proper representations of the potential elastic displacement field, we recur to the Helmholtz theorem that is commonly used to express d(x,y,z) as the sum of the gradient \(\nabla\) of a scalar potential \(\phi\)(x,y,z) and the curl \(\nabla\:\times\) of a vector potential \(\varvec{\upphi}\)(x,y,z), with components \({\upphi}_{x}\)(x,y,z), \({\upphi}_{y}\)(x,y,z) and \({\upphi}_{z}\)(x,y,z) [e.g., 32].

In the elastic regime, d(x,y,z) is a potential field if \(\nabla\:\times\:\varvec{\upphi}\)(x,y,z) is vanishing and u(x,y,z), v(x,y,z), w(x,y,z) become functions of only \(\phi\)(x,y,z), as follows [e.g., 32]:

$$\left[\begin{array}{c}u(x,y,z)\\\:v(x,y,z)\\\:w(x,y,z)\end{array}\right]=\left[\begin{array}{c}\frac{\partial\:}{\partial\:x}\phi(x,y,z)\\\:\frac{\partial\:}{\partial\:y}\phi(x,y,z)\\\:\frac{\partial\:}{\partial\:z}\phi(x,y,z)\end{array}\right]$$
(1)

From Eq. (1), it follows that u(x,y,z), v(x,y,z), and w(x,y,z) can be integrated to calculate \(\phi\)(x,y,z) over the spatial boundaries in which they are defined, here respectively indicated as xmin and xmax, ymin and ymax, zmin and zmax; for the sake of simplicity, we use the symbolism \({\phi}_{u}(x,y,z)\), \({\phi}_{v}(x,y,z)\) and \({\phi}_{w}(x,y,z)\), instead of \(\phi\)(x,y,z), as follows:

$$\begin{array}{c}{\phi}_{u}(x,y,z)={\int\:}_{{x}_{min}}^{{x}_{max}}u(x,y,z)dx\\\:{\phi}_{v}(x,y,z)={\int\:}_{{y}_{min}}^{{y}_{max}}v(x,y,z)dy\\\:{\phi}_{w}(x,y,z)={\int\:}_{{z}_{min}}^{{z}_{max}}w(x,y,z)dz\end{array}$$
(2)

where, from Eq. (1):

$$\phi\left(x,y,z\right)={\phi}_{u}\left(x,y,z\right)={\phi}_{v}\left(x,y,z\right)={\phi}_{w}(x,y,z)$$
(3)

By substituting Eq. (2) into Eq. (1), any field component can be expressed in terms of one of the others, as follows:

$$\left[\begin{array}{c}u(x,y,z)\\\:v(x,y,z)\\\:w(x,y,z)\end{array}\right]=\left[\begin{array}{c}\frac{\partial\:}{\partial\:x}\left({\int\:}_{{y}_{min}}^{{y}_{max}}v(x,y,z)dy\right)\\\:\frac{\partial\:}{\partial\:y}\left({\int\:}_{{x}_{min}}^{{x}_{max}}u(x,y,z)dx\right)\\\:\frac{\partial\:}{\partial\:z}\left({\int\:}_{{x}_{min}}^{{x}_{max}}u(x,y,z)dx\right)\end{array}\right]=\left[\begin{array}{c}\frac{\partial\:}{\partial\:x}\left({\int\:}_{{z}_{min}}^{{z}_{max}}w(x,y,z)dz\right)\\\:\frac{\partial\:}{\partial\:y}\left({\int\:}_{{z}_{min}}^{{z}_{max}}w(x,y,z)dz\right)\\\:\frac{\partial\:}{\partial\:z}\left({\int\:}_{{y}_{min}}^{{y}_{max}}v(x,y,z)dy\right)\end{array}\right]$$
(4)

Specifically, any of u(x,y,z), v(x,y,z), w(x,y,z) can be derived from the others after using directional integration and differentiation.

d(x,y,z) is not a potential field if \(\nabla\:\times\:\varvec{\upphi}\)(x,y,z) ≠ 0 and Eq. (1), Eq. (2) and Eq. (3) no longer apply. Nevertheless, two of the three components of \(\nabla\:\times\:\varvec{\upphi}\)(x,y,z) may vanish or also be isotropic. It follows that at least two of the three components of d(x,y,z) can be represented as the gradient of the same function. By considering the case of the horizontal directions under these conditions, Eq. (3) can be rewritten expressing the following condition:

$${\phi}_{u}\left(x,y,z\right)={\phi}_{v}\left(x,y,z\right)\ne\:{\phi}_{w}(x,y,z)$$
(5)

and, in turn, Eq. (4) becomes:

$$\left[\begin{array}{c}u(x,y,z)\\\:v(x,y,z)\\\:w(x,y,z)\end{array}\right]\begin{array}{c}=\\\:=\\\:\ne\:\end{array}\left[\begin{array}{c}\frac{\partial\:}{\partial\:x}\left({\int\:}_{{y}_{min}}^{{y}_{max}}v(x,y,z)dy\right)\\\:\frac{\partial\:}{\partial\:y}\left({\int\:}_{{x}_{min}}^{{x}_{max}}u(x,y,z)dx\right)\\\:\frac{\partial\:}{\partial\:z}\left({\int\:}_{{x}_{min}}^{{x}_{max}}u(x,y,z)dx\right)\:\end{array}\right]\begin{array}{c}\ne\:\\\:\ne\:\\\:=\end{array}\left[\begin{array}{c}\frac{\partial\:}{\partial\:x}\left({\int\:}_{{z}_{min}}^{{z}_{max}}w(x,y,z)dz\right)\\\:\frac{\partial\:}{\partial\:y}\left({\int\:}_{{z}_{min}}^{{z}_{max}}w(x,y,z)dz\right)\\\:\frac{\partial\:}{\partial\:z}\left({\int\:}_{{y}_{min}}^{{y}_{max}}v(x,y,z)dy\right)\:\end{array}\right]$$
(6)

,

where only u(x,y,z) and v(x,y,z) can be derived from each other using directional integration and differentiation procedures, while w(x,y,z) can neither be correctly derived from u(x,y,z) and v(x,y,z) nor be used to compute them.

Synthetic tests on simulated displacement data

Here, we show the validity and the usefulness of the proposed relations by applying them to simulated d. We use COMSOL Multiphysics software to employ a Finite Element (FE) environment in modeling sources and half-spaces with variable features. The simulated datasets have been corrupted with zero-mean normally distributed random noise terms with a standard deviation of 0.5 cm. We then apply integration and differentiation operations to the simulated field components to evaluate \({\phi}_{v}\), \({\phi}_{u}\), \({\phi}_{w}\) and, in turn, v from \({\phi}_{u}\) and \({\phi}_{w}\). Directional integration and vertical differentiation are performed in the wavenumber domain through the discrete Fourier Transform [e.g.,44]; horizontal derivatives are instead computed using finite-difference relations in the space domain. We specify that the vertical integration procedure is more robust to noise than the horizontal one; accordingly, a low-pass filter is used when the integration along the horizontal direction is carried out. Further details on integration and differentiation procedures are reported in the S1 section of the Supplementary Materials.

We start from the reference case where d satisfies Eq. (1). In volcanic framework, Eq. (1) is verified when d is produced by hydrostatically pressurized source embedded in homogeneous media [e.g., 32,33,34,35]; further, it has been shown that the potential theory is still applicable in the case of layered half-spaces17. We therefore simulate a deformation event produced by a 5 MPa hydrostatically over-pressurized (Px=Py=Pz) sill-like source with irregular geometry embedded in a layered medium with ranging Young’s Modulus (E) from 7 to 20 GPa and constant Poisson’s ratio of 0.25 (-); the modeled source extends from a depth of -5000 to -3000 m with a length of approximately 15,000 m and a width ranging from 2000 to 8000 m. The model setup is represented in Fig. 1a. We apply the integration operators to the simulated u, v, w, as described in Eq. (2), in order to compute the related functions \({\phi}_{u}\) (Fig. 1b), \({\phi}_{v}\) (Fig. 1c) and \({\phi}_{w}\) (Fig. 1d). In particular, we observe that the condition expressed by Eq. (3) is verified since \({\phi}_{u}\approx\:{\phi}_{v}\approx\:{\phi}_{w}\) (Fig. 1b, c,d). We therefore differentiate \({\phi}_{u}\) (Fig. 1b) and \({\phi}_{w}\) (Fig. 1d) along the y-direction and compare the retrieved N-S displacement from both \({\phi}_{u}\) (Fig. 1e) and \({\phi}_{w}\) (Fig. 1f) with the expected (i.e., simulated) v (Fig. 1g): the simulated and the calculated components match very well. Further details about the simulated components and the residuals analysis are reported in Figure S1 of the supplementary materials.

Fig. 1
figure 1

Synthetic tests: hydrostatic source case. (a) Model setting of the simulated scenario where the red colored body represents the hydrostatically over-pressurized source embedded in a layered half-space with variable density (ρ) and Young’s Modulus (E) parameters. Maps of the computed functions (b) \(\phi\)u, (c) \(\phi\)v and (d) \(\phi\)w. (e) N-S displacement computed from \(\phi\)u and (f) computed from \(\phi\)w; (g) Expected v.

Full size image

The second simulation is representative of a scenario described by Eq. (5). Similarly to the previously described argumentations about the potential displacement fields, we suppose that any case of sources embedded in half-spaces not causing strong anisotropies along the x- and y-directions belongs to the particular case for the application of the potential theory to the displacement field. We test this hypothesis by simulating a multi-source case where a negative pressure variation within a sill-like source simultaneously occurs with a positive one within an above pipe-like body, both embedded in a layered medium with ranging Young’s Modulus (E) from 7 to 20 GPa and constant Poisson’s ratio of 0.25 (-). In particular, the sill-like body is simulated by applying a negative 40 MPa pressure variation to the top and bottom of a 4000 × 2500 × 200 m3 rectangular prism, with − 3100 m depth to the center and a 30° strike angle with respect to the x-direction; the pipe-like source is instead simulated by applying a positive 55 MPa pressure variation to the lateral face of a vertical cylinder with radius and height respectively equal to 250 and 2000 m, and depth to the center of -2000 m. The model setup is represented in Fig. 2a. We apply the integration operators to the simulated u, v, w and calculate \({\phi}_{u}\) (Fig. 2b), \({\phi}_{v}\) (Fig. 2c), \({\phi}_{w}\) (Fig. 2d), noting that Eq. (5) is now verified since \({\phi}_{u}\approx\:{\phi}_{v}\ne\:{\phi}_{w}\) (Fig. 2b, c,d). Finally, we differentiate \({\phi}_{u}\) (Fig. 2b) and \({\phi}_{w}\) (Fig. 2d) along the y-direction and compare the N-S displacement retrieved from both \({\phi}_{u}\) (Fig. 2e) and \({\phi}_{w}\) (Fig. 2f) with v (Fig. 2g): the calculated component from \({\phi}_{u}\) well matches with the simulated one, while the use of \({\phi}_{w}\) leads to a different N-S displacement evaluation. Further details about the simulated components and the residuals analysis are reported in Figure S2 of the supplementary materials.

Fig. 2
figure 2

Synthetic tests: not hydrostatic source case. (a) Model setting of the simulated scenario where the red and orange colored bodies represent displacement sources with positive and negative pressure variations, respectively, embedded in a layered half-space with variable density (ρ) and Young’s Modulus (E) parameters. Maps of the computed functions (b) \(\phi\)u, (c) \(\phi\)v and (d) \(\phi\)w. (e) N-S displacement computed from \(\phi\)u and (f) computed from \(\phi\)w; (g) Expected v.

Full size image

Workflow diagram of the developed methodology

According to the described theoretical argumentations and performed simulations, we propose a methodological approach for the inferring of the 3-D ground displacement through three steps, as described in the workflow diagram in Fig. 3. The first phase is representative of the InSAR processing that leads to ground displacement measurements: the precise evaluation of the E-W (i.e., u) and the U-D (i.e., w) components of the ground displacement field is performed using multi-viewing (ascending/descending) SAR data. The retrieved InSAR-driven 2-D ground displacement field is then used in the second phase: \({\phi}_{u}\) and \({\phi}_{w}\) are computed and analyzed to investigate two scenarios: scenario (I) can be representative of the general solution of the potential displacement field since \({\phi}_{u}\) and \({\phi}_{w}\) are equal, and Eq. (3) becomes to be satisfied if \({\phi}_{v}\) is assumed to be also equal to them; scenario (II) can be instead representative of the particular case of the application of the potential theory to the displacement field when \({\phi}_{u}\) and \({\phi}_{w}\) are not equal, and Eq. (5) becomes true if \({\phi}_{u}\) and \({\phi}_{v}\) are assumed to be equal. Finally, the third phase returns as output the displacement along the N-S direction, i.e., v. For (I), v is computed from both \({\phi}_{u}\) and \({\phi}_{w}\), while for (II) v is calculated only from \({\phi}_{u}\), completing the 3-D ground displacement dataset.

Fig. 3
figure 3

Workflow diagram. Flowchart of the proposed approach to retrieve the 3-D ground displacement velocity maps and time series by exploiting InSAR measurements and potential theory.

Full size image

Error budget analysis on simulated displacement data

In this section, the results of some preliminary experiments carried out on simulated datasets are shown, addressing the precision of the obtained 3-D ground displacement products. For this purpose, starting from the simulated \({\phi}_{u}\) and \({\phi}_{v}\), we computed the relevant N-S components of the ground deformation field and compared their modeled values with those retrieved by applying the proposed algorithm. Specifically, we explored the following two cases: the first one in which the input modelled 2-D data (i.e., the Up-Down and East-West ground displacements) are corrupted by additive (Gaussian) noise with known statistics; in the second one the potential field assumption of hydrostatic pressure is partially relaxed to evaluate the robustness of the procedure.

Some remarks regarding the error budget of ground displacement products are preliminarily addressed. Let us refer to the acquisition geometry depicted in Fig. 4 and assume that the ascending and descending illumination geometries have approximately the same side-looking angle (i.e., \(\theta\:={\theta\:}_{asc}={\theta\:}_{Desc})\). For simplicity, let us assume the availability of N ascending and N descending SAR images independently collected in homologous dates (i.e., the time vector is \(\varvec{t}={\left[{t}_{1},\:{t}_{2},\dots\:,{t}_{N},\right]}^{T}\)). Under these simplified assumptions, we can calculate the up/down and the east/west ground deformation time series as follows:

$$\text{E}\text{a}\text{s}\text{t}\left({t}_{i}\right)=\frac{Asc.\left({t}_{i}\right)+Desc.\left({t}_{i}\right)}{2\text{c}\text{o}\text{s}\left(\theta\:\right)}\:\:i=\text{1,2},\dots\:,\:N.$$
(7)
$$\text{U}\text{p}\left({t}_{i}\right)=\frac{Asc.\left({t}_{i}\right)-Desc.\left({t}_{i}\right)}{2\text{s}\text{i}\text{n}\left(\theta\:\right)}\:\:i=\text{1,2},\dots\:,\:N.$$
(8)

where \(Asc.\) and \(Desc.\) are the LOS-projected ground displacement time series retrieved by applying the SBAS algorithm38 to the ascending and descending datasets, separately. In particular, the \(Asc.\) and the \(Desc.\) deformation time series vectors are obtained by solving the following systems of linear equations:

$$B\cdot\:Asc.=\frac{\lambda\:}{4\pi\:}\varvec{\Delta\:}{\varvec{\varphi}}_{\varvec{A}\varvec{s}\varvec{c}.}$$
(9)
$$B\cdot\:Desc.=\frac{\lambda\:}{4\pi\:}\varvec{\Delta\:}{\varvec{\varphi}}_{\varvec{D}\varvec{e}\varvec{s}\varvec{c}.}$$
(10)

where \(B\) is the design matrix of the adopted linear transformation and \(\varvec{\Delta\:}{\varvec{\varphi}}_{\varvec{A}\varvec{s}\varvec{c}.}\)and \(\varvec{\Delta\:}{\varvec{\varphi}}_{\varvec{D}\varvec{e}\varvec{s}\varvec{c}.}\) are the vectors of the M unwrapped phases related to the ascending and descending datasets, respectively.

Fig. 4
figure 4

Illumination geometries of the ascending and descending SAR tracks.

Full size image

Considering the basic rules of error propagation, the covariance matrices of the generated \(Asc.\) and \(Desc.\) time series can be expressed as follows:

$${\varvec{C}}_{\mathbf{A}\mathbf{s}\mathbf{c}.}={{\left(\frac{\lambda\:}{4\pi\:}\right)}^{2}(\varvec{B}^{\dag}}{\cdot\:\varvec{C}}_{\varvec{\Delta\:}{\varvec{\varphi}}_{\varvec{A}\varvec{s}\varvec{c}.}}\cdot\:{\varvec{B}}^{{\dag}^{\varvec{T}}})$$
(11)
$${\varvec{C}}_{\mathbf{De}\mathbf{s}\mathbf{c}.}={{\left(\frac{\lambda\:}{4\pi\:}\right)}^{2}(\varvec{B}^{\dag}}{\cdot\:\varvec{C}}_{\varvec{\Delta\:}{\varvec{\varphi}}_{\varvec{De}\varvec{s}\varvec{c}.}}\cdot\:{\varvec{B}}^{{\dag}^{\varvec{T}}})$$
(12)

where \({\varvec{C}}_{\varvec{\Delta\:}\varvec{\varphi}}\) is the covariance matrix of the unwrapped phases, and \({\varvec{B}}^{†}\) is the generalized inverse matrix of \(\varvec{B}\). Furthermore, under the simplified hypothesis that the unwrapped phases are uncorrelated, \({\varvec{C}}_{\varvec{\Delta\:}{\varvec{\varphi}}_{\varvec{A}\varvec{s}\varvec{c}.}}={\sigma}_{Asc.}^{2}{\varvec{I}}^{\varvec{M}\times\:\varvec{M}}\) and \({\varvec{C}}_{\varvec{\Delta\:}{\varvec{\varphi}}_{\varvec{D}\varvec{e}\varvec{s}\varvec{c}.}}={\sigma}_{Desc.}^{2}{\varvec{I}}^{\varvec{M}\times\:\varvec{M}}\), where \({\varvec{I}}^{\varvec{M}\times\:\varvec{M}}\) is the identity matrix. Considering that the Asc. and Desc. time series are uncorrelated, we can eventually compute the general expression for the covariance matrices of the up/down and east/west deformation time series as follows:

$${\mathbf{C}}_{\text{U}-\text{D}}=\frac{{\varvec{C}}_{\mathbf{A}\mathbf{s}\mathbf{c}.}+{\varvec{C}}_{\mathbf{D}\mathbf{e}\mathbf{s}\mathbf{c}.}}{4{\text{cos}}^{2}\left(\theta\:\right)}$$
(13)
$${\varvec{C}}_{\text{E}-\text{W}}=\frac{{\varvec{C}}_{\mathbf{A}\mathbf{s}\mathbf{c}.}+{\varvec{C}}_{\mathbf{D}\mathbf{e}\mathbf{s}\mathbf{c}.}}{4{\text{sin}}^{2}\left(\theta\:\right)}$$
(14)

For the following analyses, Eqs. (13)-(14) have been simplified, assuming that the two covariance matrices or the ascending and descending datasets are the same and diagonal \({\varvec{C}}_{\mathbf{A}\mathbf{s}\mathbf{c}.}={\varvec{C}}_{\mathbf{D}\mathbf{e}\mathbf{s}\mathbf{c}.}={\sigma}^{2}{\varvec{I}}^{\varvec{N}\times\:\varvec{N}}\). Accordingly:

$${\mathbf{C}}_{\text{U}-\text{D}}=\frac{{\varvec{I}}^{\varvec{N}\times\:\varvec{N}}\:}{2{\text{cos}}^{2}\left(\theta\:\right)}\:{\sigma}^{2}$$
(15)
$${\varvec{C}}_{\text{E}-\text{W}}=\frac{{\varvec{I}}^{\varvec{N}\times\:\varvec{N}}}{2{\text{sin}}^{2}\left(\theta\:\right)}{\sigma}^{2}$$
(16)

In the presented test, we have corrupted the displacement pattern related to the above modeled hydrostatically over-pressurized (Px=Py=Pz) sill-like source, embedded in a layered half-space (Fig. 1a) with a zero-mean 2-D Gaussian noise with covariance matrixes given by Eqs. (15)-(16), considering values of s ranging from 1 mm to 7 mm (Fig. 5a, b) and fixing the value of the incidence angle to 33°.

Fig. 5
figure 5

Error budget analysis: data accuracy. Box and whisker plots of the differences (absolute value) between the simulated (v) and computed N-S displacement from \({\phi}_{u}\) (a) and \({\phi}_{w}\) (b), for different values of σ. Blue boxes delineate the distribution from the 25th to the 75th percentile, while the red lines are the median values; red crosses represent the outliers.

Full size image

The second experiment concerns the evaluation of the error related to the assumption of the potential field theory. Specifically, we have considered the same model earlier described and simulated different cases by fixing the pressure component along x- and y-axes equal to Px=Py=15 MPa, representing the overpressure exerted to the source lateral faces, while the overpressure to the top and bottom ones has decreased from 15 to 7.5 MPa (Pz) (Fig. 6a, b). This reproduces a potential theory condition equal to \({\phi}_{w}\ne\:{\phi}_{u}\approx\:{\phi}_{v}\) (i.e., Scenario II of workflow in Fig. 3). Additionally, we have also considered the case of no longer validity of the potential theory with \({\phi}_{v}\ne\:{\phi}_{u}\approx\:{\phi}_{w}\). Here, we have fixed to Py=15 MPa the N-S oriented forces acting on the faces, while the others vary from 15 to 7.5 MPa (Px=Pz) (Fig. 6c, d).

Fig. 6
figure 6

Error budget analysis: potential field assumption. Box and whisker plots of the differences (absolute value) between the expected N-S displacement (v) and the computed from (a) \({\phi}_{u}\) and (b) \({\phi}_{w}\) versus Pz/Px = Py ratio, when the potential theory no longer verified. (c, d) Expected N-S displacement (v) and the computed from \({\phi}_{u}\) and \({\phi}_{w}\), respectively, versus Pz=Px/Py ratio. Blue boxes delineate the distribution from the 25th to the 75th percentile, while the red lines are the median values; red crosses represent the outliers.

Full size image

Results

In this section, we use the proposed method to retrieve the 3-D ground displacement field for the real case of the Sierra Negra volcano (Galapagos Islands, Ecuador). We start from the processing of 1-D InSAR data to estimate the 2-D ground displacements time-series and, finally, through the analysis of the related potential functions, infer the N-S displacements (3-D) during selected time intervals.

Sierra Negra volcano

The Sierra Negra volcano belongs to the archipelago of the Galapagos Islands (Ecuador), which constitutes the fastest deforming terrestrial volcanic province with a large present-day magma supply45. The deformation phenomena of this archipelago derive from an almost continuous magma influx through a mantle plume, producing effusive eruptions46. Sierra Negra is one of the six volcanoes occupying Isabela Island, and it is a 60 × 40 km2 basaltic shield volcano with a 7 × 9 km2 large caldera floor47. Two eruptions occurred in the past twenty years, characterized by 0.1–0.2 km3 basaltic lava flows46,48, in October 2005 and June 2018. Earthquakes (Mw > 5) and sulfur dioxide emissions have been detected associated with these events46,49.

Sierra Negra volcano hosts a continuous GNSS network composed of ten permanent stations49, of which eight fall within the caldera boundaries and are named GV01, GV03, GV04, GV05, GV06, GV07, GV08, GV09. This network was initiated in 2002 and then upgraded in 2010 to all dual-frequency receivers with a temporal sampling of 30s49. Deformation patterns, such as pre-eruptive uplift, co-eruptive subsidence, and post-eruptive uplift, have accompanied these volcanic phases49, indicating a magma migration from the mantle to shallow chamber, where it is stored, then reaching the surface with the eruptions48,49,50,51,52,53.

Geodetic data modelling of InSAR and GNSS measurements suggested that the deformation mechanisms is well represented by a ≈ 2 km shallow horizontal sill [e.g.,47,49,51] characterized by the same location and depth for at least 30 years. This finding has also been supported by gravity and tomographic models50,54.

InSAR data: february 2017 – june 2018 pre-eruptive displacement

The pre-eruptive ground displacement episodes at the Galapagos Islands have been analyzed by processing two independent sets of Sentinel-1 A SAR images collected over the investigated region and spanning the whole period between February 23, 2017, and December 28, 2018. Specifically, 65 Sentinel-1 ascending (path 106) and 65 descending (path 128) SAR images were used. The lists of available SAR data are indicated in Tables S1 and S2 of the Supplementary Material for the ascending and descending data tracks, respectively.

For every dataset, the selected images were independently co-registered to a relevant, reference common SAR image and then used to generate two sequences of small baseline (SB) SAR interferograms39,41,55. Constraints were imposed on the maximum temporal and perpendicular baseline of 100 days and 100 m, respectively, for the selected SB InSAR data pairs. For the interferogram generation, precise orbital information and the one-arc-sec SRTM DEM of the area were used56. Then, the ascending/descending (1-D) LOS projected ground displacement time-series were recovered by applying the Small Baseline Subset (SBAS) technique38,57. Finally, the 2-D (E-W and U-D) ground displacement time-series were recovered over the set of coherent points common to both ascending and descending data tracks by applying the combination methods described in Pepe et al.25 and Pepe and Calò58.

We show the retrieved mean E-W (Fig. 7a) and U-D (Fig. 7b) velocity maps related to the analyzed time interval. We also compare the retrieved E-W (red dots in Fig. 8a) and U-D (red dots in Fig. 8b) displacement time-series with the corresponding GNSS measurements (blue triangles in Fig. 8). The selected GNSS datasets are related to the stations within the caldera, i.e., GV03, GV04, GV05, GV06, GV07, GV08, GV09, whose locations are reported in Fig. 7 (black triangles); the GV01 station is used as a reference point and, therefore, is not shown in Fig. 8; also, the Root Mean Square Error (RMSE) values (indicated as \(\xi\) in Fig. 8) between InSAR and GNSS are reported. We specify that the comparison is performed by selecting the highest coherence InSAR pixel among the nearest ones to the considered GNSS station within a 200 × 200 m2 box.

GNSS datasets have been freely downloaded from the UNAVCO archive (https://www.unavco.org/data/data.html) for the analyzed period [e.g., 59] and then interpolated at each considered SAR epoch.

Fig. 7
figure 7

InSAR measurements at Sierra Negra volcano: ground displacement velocity maps. Mean (a) E-W and (b) U-D displacement velocity maps computed for the Southern region of Isla Isabela during the analyzed February 2017 – June 2018 time interval. The black box represents the zoomed area of the Sierra Negra caldera. Black triangles and white abbreviations indicate the GNSS stations’ locations and names, respectively. White squares indicate the spatial reference pixel. Coordinates are expressed in the UTM reference system (UTM zone: 15 S).

Full size image
Fig. 8
figure 8

InSAR measurements and GNSS data at Sierra Negra volcano: time-series analysis. Comparison between InSAR (red dots) and GNSS (blue triangles) time-series of the (a) E-W and (b) U-D components of the ground displacement field referred to the analyzed February 2017 – June 2018 time interval. ξ indicates the related RMSE. The GV01 GNSS station represents the reference point.

Full size image

Analysis of potentials

We have applied the integral procedures to E-W and U-D deformations at each epoch to compute the functions \({\phi}_{u}\) and \({\phi}_{w}\), respectively, and analyze their time variations. The deformation components have been gridded here using a natural neighbours interpolator with a 100 m sampling step along both directions. The grids were then extended around Isabela Island using an extrapolation method based on the smooth extension of order 0 to circumvent aliasing from circular convolution60. Finally, a low-pass filter was applied to increase the SNR for the E-W case.

The computed \({\phi}_{u}\) and \({\phi}_{w}\) are here shown in terms of statistically representative indexes for each considered InSAR time epoch, such as maxima (Fig. 9a), mean values (Fig. 9b) and variances (Fig. 9c). We observe that \({\phi}_{u}\) (orange dots in Fig. 9) and \({\phi}_{w}\) have comparable values (green dots in Fig. 9) during February 2017 – September 2017, while they diverge from September 2017 to before the eruption (June 2018). Therefore, we point out that:

  • \({\phi}_{u}\approx\:{\phi}_{w}\) condition occurs during February 2017 – September 2017 time interval; this agrees with scenario (I) that assumes \({{\phi}_{v}\approx\:\phi}_{u}\) and \({{\phi}_{v}\approx\:\phi}_{w}\); we therefore proceed by evaluating the N-S displacement from the U-D component (or the E-W one);

  • \({\phi}_{u}\approx\:{\phi}_{w}\) condition is verified for the next September 2017 – June 2018 time interval; this is instead in accordance with scenario (II) assuming \({\phi}_{v}\approx\:{\phi}_{u}\); we derive the N-S component from the E-W one.

Fig. 9
figure 9

Temporal analysis of potentials at Sierra Negra volcano. February 2017 – June 2018 time-series of the (a) maxima, (b) mean values and (c) variance of the potentials computed from the E-W (\(\phi\)u, orange dots) and U-D (\(\phi\)w, green dots) displacements from InSAR measurements. The black dashed lines mark the two identified time intervals.

Full size image

N-S displacement (i): february 2017 – september 2017 time interval

The results of differentiation procedures applied to both potentials are reported in Fig. 10. In particular, we show the N-S time-series retrieved from \({\phi}_{u}\) and \({\phi}_{w}\) (red dots in Fig. 10a and b, respectively) with the corresponding GNSS dataset (blue triangles in Fig. 10a, b). RMSE values for this comparison are reported (indicated as \(\xi\) in Fig. 10).

In Fig. 10c, we also show the mean N-S ground displacement velocity map retrieved from \({\phi}_{w}\), with values evaluated only at the coherent pixels of InSAR processing. This describes a displacement pattern affecting the entire caldera with a ~ 40 cm/yr maximum velocity towards the N and S directions.

Fig. 10
figure 10

N-S ground displacements at Sierra Negra volcano: February 2017 – September 2017 time interval. Comparison between the N-S displacement components (red dots), retrieved from (a) \(\phi\)u and (b) \(\phi\)w, and GNSS (blue triangles) time-series. ξ indicates the related RMSE. (c) Mean N-S displacement velocity map calculated from \(\phi\)w. The GV01 GNSS station represents the reference point. Coordinates are expressed in the UTM reference system (UTM zone: 15 S).

Full size image

N-S displacement (ii): september 2017 – june 2018 time interval

Here, we compute the N-S displacement for the second time interval by only applying the differentiation procedures to \({\phi}_{u}\). Also in this case, the retrieved time-series (red dots in Fig. 11a) are shown with the GNSS ones (blue triangles in Fig. 11a). RMSE are indicated as \(\xi\) in Fig. 11. Finally, the mean N-S velocity map (Fig. 11b) shows a variation of the displacement pattern within the caldera floor with a maximum velocity of ~ 60–70 cm/yr towards the N and S directions. We report also the results related to the use of \({\phi}_{w}\) in Figure S3 of the supplementary information.

Fig. 11
figure 11

N-S ground displacements at Sierra Negra volcano: September 2017 – June 2018 time interval. (a) Comparison between the retrieved N-S displacement component from \(\phi\)u (red dots) and GNSS (blue triangles) time-series. ξ indicates the related RMSE. (b) Mean N-S displacement velocity map calculated from \(\phi\)u. The GV01 GNSS station represents the reference point. Coordinates are expressed in the UTM reference system (UTM zone: 15 S).

Full size image

Discussion

The proposed method allows retrieving the 3-D ground displacement field in volcanic environments: conventional differential interferometric SAR techniques are applied to ascending/descending data tracks to estimate the E-W and the U-D ground displacement components, with uncertainties in the order of a few millimeters per year (for displacement velocity). Instead, the N-S ground displacement components are derived, with enhanced precision than that attainable using exclusively SAR/InSAR technologies, exploiting the fundamentals of potential field theory. The integration of InSAR and potential field theory looks innovative in this context.

To discuss the applicability of the proposed methodology, we use the expressions of the potential field of the elastic displacement, Eq. (1), which describe the E-W, N-S and U-D components as the gradient of the same scalar function; we then treat the particular case that only the E-W and N-S components can be expressed using the same function, according to the representative condition expressed by Eq. (5). This is fundamental for evaluating one or more components, in our case, the N-S one, from the other components using integration and differentiation procedures. For volcanic unrests, Eq. (1) is undoubtedly satisfied when sources with hydrostatic pressure change are embedded in layered elastic media; Eq. (5) is instead verified for describing with a good approximation scenario without strong anisotropies along the horizontal directions, as horizontal sills emplacement or pressurized vertical pipes embedded in layered half-spaces.

The application of potential theory to simulated data confirmed our argumentations for volcanic contexts. In the first test (Fig. 1), i.e., when hydrostatically pressurized sources with any geometry are embedded in elastic layered media, the N-S component was well derived from both E-W and U-D displacements with discrepancies of a few millimeters with respect to the expected one (Figure S1c, d for further details). This was instead true only for the E-W component for the second test (Fig. 2), where a magma migration phenomena from a sill- to pipe-like sources in a layered medium is simulated; also in this case, the residual analysis showed that the N-S component is evaluated with discrepancies of a few millimeters with respect to the expected ones (Figure S2c for further details). These tests were performed under controlled conditions, we can associate these discrepancies with using Fast Fourier Transform algorithms for performing integration and differentiation procedures or the instability of the horizontal integration versus noise. Both tests proved positive since the maximum intensity of the retrieved discrepancies is lower than the noise used to corrupt the simulated fields.

To clarify which component has to be used for a more reliable estimate of the N-S displacement, we proposed a flowchart involving two analysis scenarios, after estimating the 2-D ground displacement field from InSAR measurements (Fig. 3). In the first scenario (I), it is verified that the same potential characterizes the E-W and U-D components, and the N-S one can be evaluated from both, assuming that it also shares the same potential, in accordance with Eq. (3). We specify that, in areas where the incidence angle on the ground of the radar pulses is significantly smaller than 45°, the use of the U-D ground displacement is preferred, since it is more precise than the E-W component. In the second scenario (II), E-W and U-D components have different potentials, and the N-S deformation is estimated from the E-W one, assuming that only the latter is the gradient of the same potential, as implied by Eq. (5).

Although the analysis of potential functions may provide indications as to the reliability of the retrieved N-S displacement from the E-W and U-D components, the use of the potential theory does not allow an unambiguous quantification of the error related to the inferred displacements without recurring to independent information, such as GNSS data. Nevertheless, the proposed preliminary assessment of the expected errors during the displacement retrieval attempted to clarify the potentiality and limitation of this approach. The first test showed an increase in the error in retrieving the N-S displacement as the accuracy of input InSAR-derived E-W (Fig. 5a) and U-D (Fig. 5b) components worsened. The test confirmed that the more precise these measurements, the fewer the errors associated with the use of the potential theory procedure. For larger σ value, i.e., from 0.5 to 0.7, the use of the E-W displacements involves larger errors in the N-S retrieval, which can also reach about 2 cm in maximum discrepancy (Fig. 5a) than the 1 cm related to the use of the U-D one (Fig. 5b). Considering a real case, these errors are acceptable for such task. The second test, related to the validity of potential theory assumptions, pointed out that the E-W displacement allows the retrieval of the N-S component with lower error than 0.2 cm as the anisotropy of pressure between the x-y and z-directions increases (Fig. 6a). This error can instead reach maxima of about 3 cm considering the U-D one (Fig. 6b), whose use is therefore not recommended for these scenarios. Finally, the third case highlighted that errors in retrieving the N-S displacements from both the E-W (Fig. 6c) and U-D (Fig. 6d) components can also exceed 6 cm in case of stresses anisotropy between the x-z and y-directions.

The use of potential theory is therefore not recommended to analyze the cases of strong variations of the pressure and so the stress tensor since they can lead to errors that cannot be unambiguously quantified. This can be the case of an intruding dike with a preferential opening along the N-S directions or of any source embedded in a strongly anisotropic half-space in its immediate surroundings. It follows that the use of potential theory is also discouraged to study phenomena governed by strong shear stresses, as earthquakes or landslides.

For the real case of the Sierra Negra volcano, we propaedeutically addressed the quality of the processed time-series through a comparison with GNSS measurements, showing a maximum discrepancy between them, in terms of RMSE, of about 1 cm (Fig. 8). Independent GNSS data proved helpful in evaluating the performance of the potential theory and associating a rough error to the estimates of N-S ground displacement components. Specifically, we got a mean RMSE values of ~ 1.3 cm and ~ 1.75 cm (maximum values of 2 and 3.5 cm, respectively) in the case of early February 2017 – September 2017 pre-eruptive phase (Fig. 10) and in the eruption proximity during September 2017 and June 2018 (Fig. 11). With this kind of validation, our approach becomes a reliable tool that can provide spatially denser information with respect to the GNSS dataset, extending the knowledge of N-S displacement from less than ten GNSS stations to tens of thousands of measurement points.

The analysis of the potential functions (phase 2 of the proposed workflow) at Sierra Negra volcano allowed the identification of two different intervals of the analyzed time-series, where shape variations of the active displacement source could have occurred from February 2017 – September 2017 to September 2017 – June 2018 since the two scenarios proposed in the workflow implicitly involve different volcanic mechanisms. This argumentation is valid since the eruption occurred at the end of the analyzed period, i.e., June 2018, and changes and/or migration of the displacement sources are expected from the pre-eruptive to the eruptive phases. We may, therefore, argue that the first interval was characterized by a stationary magma chamber with hydrostatic pressure variations and that the magma subsequently rose during the second interval, causing the eruption on 26th June 2018. Therefore, it is interesting to note that our method could be helpful for monitoring purposes, helping understand when a variation of the displacement mechanism occurs during a volcanic unrest.

Interesting future developments are related to this work, such as: (i) the impact evaluation of using a 3-D ground deformation field in the modeling framework; (ii) the monitoring of the source-mechanism evolution before eruptive phases using the potential theory; (iii) the study of phenomena governed by the poro-elasticity laws using potential theory61. In case (i), we have seen, in the introduction section, how the use of the 3-D ground deformation dataset could constrain deformation sources during modeling procedures, helping to achieve more reliable solutions; however, the proposed strategy involves a dependence on the N-S deformation on the others, and, therefore, it is necessary to confirm the contribution of our 3-D field for modelling purposes. In case (ii), the application to the real case of Sierra Negra volcano showed how time-varying potentials can reveal changes in the deformation sources, and it would be interesting to understand the potentiality of this tool for monitoring the source shape variations, which can often prelude an eruption. Finally, in case (iii), the potential theory can also be used to study phenomena governed by the poro-elasticity law, such as deformations due to aquifer dynamics or geothermal fields; however, this case has not been explored in this work.