# Mitigating bias in inversion of InSAR data resulting from radar viewing geometries

Quentin Dumont, Valérie Cayol, Jean-Luc Froger

► To cite this version:

Quentin Dumont, Valérie Cayol, Jean-Luc Froger. Mitigating bias in inversion of InSAR data resulting from radar viewing geometries. Geophysical Journal International, 2021, 227 (1), pp.483-495. 10.1093/gji/ggab229. hal-03356622

HAL Id: hal-03356622

https://uca.hal.science/hal-03356622

Submitted on 28 Sep 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Mitigating bias in inversion of InSAR data resulting from radar viewing geometries

Quentin Dumont¹, Valérie Cayol¹, Jean-Luc Froger¹

¹ Laboratoire Magmas et Volcans, OPGC, Université Clermont-Auvergne, CNRS, IRD, 6 av. Blaise Pascal, 63170 Aubière, France

Volume 227, Issue 1, October 2021, Pages 483–495, https://doi.org/10.1093/gji/ggab229

Published 2021 June 16; Accepted 2021 June 14. Received 2021 June 04; in original form 2020 October 29

## SUMMARY

InSAR data acquired from ascending and descending orbits are often characterized by different magnitudes of the observed line-of-sight displacements, which may potentially bias inverse models. Using synthetic numerical models of dyke intrusions, we show that biased solutions are obtained when carrying out "conventional" inversions where only observation and modelling errors are taken into consideration. To mitigate the impact of the relative magnitudes of the data, we propose two methods: a covariance weighting inversion and a wrapped data inversion. These methods are compared to a conventional inversion using synthetic data generated by models of dykes of known geometry. We find that the covariance weighting method allows to retrieve an initial source geometry better than the other methods. These methods are then applied to the July 2017 eruption of Piton de la Fournaise. Using a covariance weighting inversion, the difference in fit between data sets decreases from 50% to 20% and the newly estimated source is in better agreement with the geological context.

**Key words:** Radar interferometry – Indian Ocean – Interferometry – Inverse theory – Numerical modelling – Remote sensing of volcanoes.

## 1 Introduction

InSAR measures ground surface displacement projected along a satellite's line-of-sight (LOS), which depends on the satellite heading (ascending or descending), radar direction (left or right looking) and acquisition mode. Because of this projection, the measured displacement is inherently ambiguous and the same 3-D displacement vector can be represented by completely different LOS values. When a limited number of LOS is available, the 3-D displacement field is incompletely characterized, which might affect the determination of the sources of deformation through inverse modelling (Wright et al. 2004). Acquisitions in several LOS in ascending and descending orbits make the displacement field and the characterization of the deformation sources less ambiguous (Fialko et al. 2001; Wright et al. 2004; Dawson & Tregoning 2007; Remy et al. 2014; Peltier et al. 2017).

When different types of data are used, such as GNSS and InSAR, in which the number of measurements and the uncertainties are also different, the question is raised of the data weighting required for an accurate determination of the sources of displacement. Given that the number of sample points of GNSS is less than that of InSAR, GNSS data have little weight on the inversion result when only the intrinsic data uncertainties are taken into consideration, as advocated by Tarantola (2005) and applied for instance by Árnadóttir et al. (2004). Recognizing that this approach might not be sufficient to counterbalance uneven numbers of points in data sets, several weighting strategies have been developed. Fukuda & Johnson (2008, 2010) estimated the optimal data weighting simultaneously with the model parameters. Fialko (2004) and Auriac et al. (2013) weighted the covariance proportionally to the number of points of each data sets and Smittarello et al. (2019b) gave the same likelihood to InSAR and GNSS data by normalizing the covariance matrix of each data set.

Here, following Smittarello et al. (2019b), we argue that the number of points is not the only factor that may affect the outcome of an inversion, but that the data magnitude may also have a significant influence. This issue has not been properly investigated up to now. For Piton de la Fournaise dyke intrusions, inversions of InSAR displacement often result in models providing a higher percentage of
<figure><img src="image_2.png" /><figcaption>Figure 1. Relative weight of the descending data over the total displacements recorded, $W_{Dsc}$ (eq. 1). The large squares indicate the dyke configurations tested in this study. The pink vertical bar indicates the range of dip values expected at Piton de la Fournaise volcano.</figcaption></figure>

explained data on the descending pass than on the ascending one (Froger et al. 2004; Fukushima et al. 2010; Smittarello et al. 2019b). For example, the models presented in Smittarello et al. (2019b) typically explained descending data sets at 94-97%, while ascending data sets were only explained at 57%-87%. This has generally been attributed to an incorrect mechanical framework as postulated by Letourneur et al. (2008); Got et al. (2013); Carrier et al. (2015), but, here, we postulate that it could also be attributed to a preferential convergence toward the data set with the highest average magnitude.

As LOS projections can lead to differences in data magnitude, we investigate their impact on inversions. We generated synthetic InSAR data from 3D dykes submitted to a constant overpressure by using a 3-D Mixed Boundary Element method (3D-BEM) for elastic, linear and homogeneous media (Cayol & Cornet 1997, 1998). Realistic topographies are taken into account. Dykes and the ground surface are meshed with triangular elements. Following estimations of Cayol & Cornet (1998) for Piton de la Fournaise and consistent with in-situ measurements of Heap et al. (2020) for basaltic volcanoes, the medium was assumed to have a Young's modulus of 5 GPa and a Poisson's coefficient of 0.25. Simulated dykes are rectangular, 500 m high by 800 m long, and submitted to an overpressure of 0.5 MPa. Four strikes, N30°E, N90°E, N140°E and N180°E, and several dips, vertical to 20° to the east or west were considered. Synthetic interferograms were then produced by projecting the obtained 3D displacements in ascending and descending directions using the typical LOS of the Sentinel-1 satellite at La Réunion Island (see LOS in Supporting Information Table S2).

As a proxy for the relative weight of a data set during an inversion, we calculate the ratio between the descending data $U_{Dsc}$ and the sum of the ascending and descending data $U_{Dsc+Asc}$:

$$W_{Dsc} = \frac{U_{Dsc}^T U_{Dsc}}{U_{Dsc+Asc}^T U_{Dsc+Asc}}, \quad \text{where} \quad U_{Dsc+Asc} = \begin{bmatrix} U_{Dsc} \\ U_{Asc} \end{bmatrix} \quad (1)$$

where $U_{Dsc}$ and $U_{Asc}$ are vectors containing all displacement values of the descending and ascending interferograms, respectively. $W_{Dsc}$ corresponds to the contribution of the descending interferogram to the total displacements. Values of $W_{Dsc}$ (see Fig. 1) show that, for all the N90°E striking dykes, $W_{Dsc} \approx 0.5$, indicating that descending and ascending interferograms have the same weight. For other strikes, data relative weights can vary from less than 0.05 to more than 0.95, with maximum and minimum values of $W_{Dsc}$ corresponding to dips of 60°.

At Piton de la Fournaise volcano, typical dykes in the rift zones show strikes varying from N140°E to N30°E with eastward dips between 60°E and 90°E (Froger et al. 2004; Fukushima et al. 2010) as represented in Fig. 1 by the pink bar. Only a few eruptions are fed by N90°E dykes. In this type of field configuration, unbalanced weights of the data sets may be critical when estimating dyke geometries. In a context where intrusive systems are oriented N90°E, such as for the Kilauea East Rift Zone, the influence of the magnitude ratio should be negligible.

In this paper, we investigate the bias caused by differences in the magnitude of LOS displacements on inverted sources as well as ways to compensate for them. First, we show the importance of the data magnitudes on the inversion results and propose two strategies to counteract potential biases. Their effectiveness is assessed through synthetic tests. As a test case, we apply and compare these methods for the inversion of the displacement measured during the July 2017 eruption of Piton de la Fournaise.
## 2 Methods

### 2.1 InSAR data inversions

Inverse modelling relies on observations to infer the values of the parameters characterizing a system (Tarantola 2005). We can express the displacement observations as a function of the model parameters as:

$$u_{obs} = F(m) + \epsilon \quad (2)$$

where $u_{obs}$ is the observed displacement, $F$ is the forward modelling operator, which is a non linear function of the source location and geometry, $m$ is the vector of input parameters, and $\epsilon$ is the vector of measurements and modelling errors. When inverting InSAR data, measurement errors are generally linked to atmospheric effects (Zebker et al. 1997; Doin et al. 2009) or variable time spans of the interferometric pairs, while the modelling errors come from inexact assumption about the source types, number of sources (Chen et al. 2017; Smittarello et al. 2019b), inadequate evaluation of the source interactions (Pascal et al. 2013), or an inexact mechanical framework (Letourneur et al. 2008; Montgomery-Brown et al. 2011; Got et al. 2013; Duputel et al. 2015).

Inverting the sources of displacements consists of searching for the value of $m$ that minimizes the difference between the model and the observations, within given measurements and model errors. It is defined by the cost function as:

$$\chi^2 = (u_{obs} - u_{mod})^T C_d^{-1} (u_{obs} - u_{mod}) \quad (3)$$

where $u_{obs}$ and $u_{mod}$ are respectively the vector of observed and modelled displacements, and $C_d$ is the covariance matrix of the data which makes it possible to weight data according to their variance.

As InSAR data show correlated random noise, we follow Tarantola (1987) for the covariance calculation of each data set. Fukushima et al. (2005) have shown that InSAR noise correlation decreases exponentially with distance. We follow their formulation of the covariance matrix defined as:

$$C_d(d) = \sigma^2 \times e^{-d/a} \quad (4)$$

where $d$ is the distance between two points, $\sigma^2$ is the noise variance, and $a$ is the correlation length. As shown by Tarantola (2005), the variance $\sigma^2$ results from (1) errors in the data and (2) errors of the model ($C_d = C_{interferogram} + C_{model}$). This variance is estimated through the residuals of preliminary inversions. For Piton de la Fournaise interferograms covering five eruptions, Fukushima et al. (2010) use best-fitting model residuals to determine that a variance of $5 \times 10^{-4} m^2$ and a correlation length of 850 m are representative of both data and modelling errors. These values are the one taken for this study.

In contrast to forward modelling, where a set of parameters results in a unique solution, the inverse problem is ill-posed, leading to several models and parameter values that explain the data equally well. In consequence, any a priori constraints on the models are valuable to better determine the relevant model and restrict parameter values. They can be provided by field data, or other geophysical measurements. Alternatively, continuous GNSS data may allow discriminating between equally likely models derived from static InSAR data as shown by Smittarello et al. (2019b).

The range of possible values for each parameter of $m$ is called the parameters space. Because ground displacements are a non linear function of the parameters, the inversion of the source location and geometry consists of exploring the parameter's space following an a priori probability density function, to determine the range of parameter values that explain the data within their uncertainties. In our case, we assumed the a priori probability density function to be uniform.

Inversions are conducted in two stages: a search and an appraisal stage. The search stage relies on the combination of forward model computations, with the 3D-BEM previously used to generate the synthetic displacements, and a neighborhood inversion algorithm (Sambridge 1999a; Fukushima et al. 2005), where the neighborhood is defined by Voronoi cells. Each dyke geometry is approximated by a quadrangular surface defined by nine parameters presented in the Supporting Information (Fig. S1). After calculation of the displacement for an unit overpressure of 1 MPa, the dyke overpressure is scaled to minimize the misfit, following the approach outlined by Tridon et al. (2016).

The inversion algorithm starts by randomly picking $N_1$ points in the parameter space, where $N_1$ is approximated by an exponential function of the number of parameters, $k$, $N_1 = 9.775 \times e^{0.457 \times k}$ (Tridon et al. 2016), which is determined by the number of natural neighbours as estimated by Sambridge (1998). For subsequent iteration, following Fukushima et al. (2005), 50 points are drawn in the neighbourhood of the 50 lowest misfit points. This enables the search to be explorative enough, while keeping the inversion time reasonable. Iterations stop when one of the following empiric conditions is reached: (1) the standard deviation on $\chi^2$ of the 50$^{th}$ last model is smaller than 0.05, (2) the scaled model parameter values of the 50$^{th}$ last model is smaller than 0.05, (3) 600 iterations have been reached. In this later case, we manually check actual convergence of the inversion. In all our inversions, satisfactory convergence was reached for 600 iterations or less (see the Supporting Information S2 and S10 for the convergence of all inversions).

In order to obtain numerically computable misfits, the interferograms are subsampled by a quadtree algorithm (Welstead 1999; Jónsson et al. 2002; Lohman & Simons 2005), to give a large number of points in high strain or displacement areas, and a small number of points
## 4 Q. Dumont, V. Cayol, J.-L. Froger

in low strain or displacement areas. In order to better account for the area each quadtree cell represents, the covariance of each cell can be computed from the full covariance matrix corresponding to all the interferogram pixels (Sudhaus & Jónsson 2009; Tridon et al. 2016). This type of covariance weighting was tested on our inversions but it did not change the results significantly. Thus, to minimize calculation time of the covariance matrix, we choose not to weight the subsampled data by the area they represent.

In the appraisal stage of the inversions, 1-D and 2-D posterior probability functions are estimated using the Bayesian inference following Sambridge (1999b). The model population calculated during the search stage is resampled by a Monte Carlo integration, allowing posterior probability density functions to be reconstructed. The confidence interval and the mean model are derived from the 1-D posterior probability density functions.

With this inversion method, hereafter referred to as the Conventional inversion, the data sets have an intrinsic weight which depends on (1) the variance of the data related to measurement and modelling errors, (2) the number of subsampled points of each data set, and (3) the magnitude of the measured displacement.

## 2.2 Methods for compensating for unbalanced weights of InSAR data

In order to give equal weights to each data set during an inversion, we investigated several methods: one involves adjusting the number of subsampled points for each data set; the second, hereafter referred to as the Covariance Weighting inversion method, involves weighting the data through the covariance matrix; and the last, hereafter referred to as the Wrapped Inversion method, involves removing the influence of data magnitude by carrying out inversions on wrapped interferograms (Feigl & Thurber 2009; Jiang & González 2020).

Adjusting the number of subsampled points for each data set to correct for the magnitude would require very different numbers of subsampled points between data sets. For instance, for the July 2017 eruption of Piton de la Fournaise, the ascending interferograms would need 12 times more points than the descending one to have both data sets contribute to the misfit value by the same amount. This would either lead to a very high number of data for the low magnitude interferogram, which would not be numerically manageable, or to a very low number of data for a high magnitude interferogram which would make the subsampling too limited to provide information on the spatial distribution of the deformation. Consequently, we do not explore this possibility further.

### 2.2.1 Covariance Weighting inversion method

Compensating for the unbalanced magnitudes by weighting the covariance matrix is straightforward. It involves multiplying the covariance matrix corresponding to each interferogram by a ratio indicative of the data set's relative magnitude (see part B of the Supporting Information for derivation details):

$$C_{d_i}^{\text{correct}} = C_{d_i} \times \left( \frac{N_d \cdot (\chi_0^2)_i}{(\chi_0^2)_{tot}} \right), \quad \text{with} \quad \chi_0^2 = \mathbf{u}^T C_d^{-1} \mathbf{u} \qquad (5)$$

where $C_{d_i}$ is the covariance sub-matrix of the i-th data set, $N_d$ is the number of data sets, $(\chi_0^2)_i$ and $(\chi_0^2)_{tot}$ are the intrinsic data weight (corresponding to a null displacement model) of the i-th data set and of all data sets, respectively. In terms of Bayesian inference, this correction gives the same likelihood to each data set.

### 2.2.2 Wrapped Inversion method

Another approach consists of carrying out inversions on wrapped data. This method has been proposed by Feigl & Thurber (2009); Ali & Feigl (2012) and improved by Jiang & González (2020) in order to overcome errors related to the phase-unwrapping step. Here, we use it in order to cancel the influence of data magnitudes. Indeed, working with phases, the variation of data magnitudes is inherently limited to $-\pi$ and $\pi$. Hence, we expect a better sensitivity to low magnitude displacements.

Wrapped Inversion uses the Angular Deviation cost function, modified from Feigl & Thurber (2009) by adding a covariance matrix term to weight data. It is defined as:

$$\Omega = \omega^T C_{d_{wrp}}^{-1} \omega \qquad (6)$$

where $C_{d_{wrp}}$ is the covariance matrix of the wrapped interferograms, and $\omega$ is a vector of angular deviations given at each point by:

$$\omega^i = \min \left( |\phi_{obs}^i - \phi_{mod}^i|, 2\pi - |\phi_{obs}^i - \phi_{mod}^i| \right) \qquad (7)$$

where $\phi_{obs}^i$ and $\phi_{mod}^i$ are the observed and modelled phases for each data i, respectively. When inverting the wrapped phase, pressure cannot be estimated via linear least-squares because phase is a non-linear function of pressure. Consequently, the pressure is inverted as a supplementary non-linear parameter. In order to compare this method with the previous ones, we worked with rewrapped data using the same subsampled points.
<table><caption>**Table 1.** Data set relative weights for different strikes and dips of synthetic dykes shown in Fig. 1. $W_{Dsc}$ is the relative weight of the descending data over the sum of the ascending and descending data (eq 1). %ED$_{Dsc}$, %ED$_{Asc}$ and %ED$_{tot}$ are the percentage of explained data for descending, ascending, and both interferograms, respectively, all unmasked points in the interferogram being taken into account. $D_{Mesh}$ gives the average distance between the inverted mesh and the initial mesh and its standard deviation, expressed in meters. Dykes are shown in Fig. 2.</caption><thead><tr><th>Dyke</th><th>$W_{Dsc}$</th><th>%ED$_{Dsc}$</th><th>%ED$_{Asc}$</th><th>%ED$_{tot}$</th><th>$D_{Mesh}$ (m)</th></tr></thead><tbody><tr><td>N140°E-60°E</td><td>0.93</td><td>95</td><td>80</td><td>92</td><td>94±144</td></tr><tr><td>N140°E-60°W</td><td>0.09</td><td>87</td><td>94</td><td>93</td><td>121±182</td></tr><tr><td>N90°E-60°S</td><td>0.52</td><td>95</td><td>94</td><td>94</td><td>45±41</td></tr></tbody></table>

## 2.3 Criteria for models comparisons

Because our objective is to have the inversions fit the observed InSAR data and to accurately determine the geometry of an intrusion, we base the comparison of each final model on three criteria.

The first is the RMS of the residuals. However, such absolute criterion does not give information on the goodness of the fit relative to the intrinsic magnitude of the signal but rather indicates the magnitude of residual. Indeed, a RMS of 1 cm on a signal of 40 cm, is indicative of a much better model than an RMS of 1 cm for a signal of 10 cm (regardless of any noise). Consequently, to estimate residuals relative to the data set average magnitude and to allow for the cross comparison of results between inversions having different cost criteria, we assess the goodness of fit of the best model by using the percentage of explained data (%ED), defined as:

$$ \%ED = \left( 1 - \sqrt{\frac{\left(\mathbf{U}_{ref} - \mathbf{U}_{mod}\right)^T \left(\mathbf{U}_{ref} - \mathbf{U}_{mod}\right)}{\mathbf{U}_{ref}^T \mathbf{U}_{ref}}} \right) \times 100 \quad (8) $$

where $U_{ref}$ and $U_{mod}$ are the reference and modelled displacements on the full interferograms (without any subsampling), respectively. This criterion indicates how well the model fits the data relative to a null reference model, expressed as $U_{ref}^T U_{ref}$.

To quantify the synthetic dyke mesh similarity, we follow the method proposed by Smittarello et al. (2019a) which consists in computing the average distance between meshes A and B, the distance between A and B corresponds to the average distance between each points of A and the points of B nearest to this point. With this criteria, the distance between A and B and between B and A are slightly different (see fig. 7 of Smittarello et al. (2019a)). Consequently, we take the mean of the both distances, noted as $D_{Mesh}$.

## 3 Synthetic Tests

In the following section, we first use synthetic tests and Conventional inversions to investigate the influence of unbalanced data weight on the data fit. We next use synthetic tests to determine a reliable inversion strategy. For this approach, we compare model characteristics and the data fit.

### 3.1 Biased inversions arising from unbalanced data weights

To test whether the relative magnitudes of data sets bias the inverted sources, we carry out Conventional inversions on three synthetic dykes: N140°E-60°E, N140°E-60°W and N90°E-60°S (see Fig. 2). These dykes are characterized by a high relative weight ($W_{Dsc} \approx 1$), a low relative weight ($W_{Dsc} \approx 0.1$) and a well-balanced relative weight ($W_{Dsc} \approx 0.5$), respectively (see Table 1 and Fig. 1). LOS displacements were computed, then inverted following the Conventional inversion method described in Section 2.1. Here, we only discuss how well the data are explained. Detailed tests on the way a dyke is biased will be conducted in the following section.

The inversion of a N140°E-60°E dyke, in which the weight of the descending data dominates ($W_{Dsc} \approx 1$), better explains the descending than the ascending one (94% versus 80%). On the contrary, the inversion for the N140°E-60°W dyke, which has a low weight on the descending data ($W_{Dsc} \approx 0.1$), better explains the ascending interferogram (87% versus 93%). Finally, the inversion of the N90°E-60°S dyke, with a well-balanced data set weight, equally explains both interferograms (94% each). In this latter case, the average distance between the initial and the inverted meshes (45 m) is significantly lower than for the N140°E-60°E and N140°E-60°W dykes (94 m and 121 m, respectively), indicating a better fit to the initial geometry when the data sets are balanced.

### 3.2 Balancing data set magnitudes

To test the effectiveness of the approaches proposed to balance data sets magnitudes, we compared inversions of synthetic interferograms corresponding to an initial dyke N140°E-60°E. A reference inversion is conducted using the Conventional inversion method. Then we tested the Covariance Weighting method, and the Wrapped inversion method. We also compared results with inversions conducted
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Dykes and line-of-sight (LOS) displacements considered in the synthetic tests used to test potential bias the data weight has on inversions. (a.) Dykes geometries and openings in map view. (b.) Displacements projected along descending and ascending Sentinel-1 LOS (see Supporting Information Table S2). For each plot, thick black lines indicate contours of -0.005 and 0.005 m, while dashed lines indicate a contour equal to zero. Black arrows in the bottom corner of each panels indicate LOS directions.</figcaption></figure>

using either a descending interferogram or an ascending interferogram alone (referred to as Conventional Descending and Conventional Ascending inversions, respectively). Moreover, to ensure comparability between inversions, the same subsampling points were used with 429 and 432 points for the descending and ascending interferogram, respectively (Figs. S3 to S7 in the Supporting Information). Please note that no noise was added to the synthetic interferograms to isolate the influence of the intrinsic magnitudes. Because of the highly non-linear geometry description by parameter values, a given confidence interval on one parameter does not have necessarily the same geometrical implication when associated to a different set of other parameters. We have therefore not discussed directly the values spanned by the 95% confidence intervals (Supporting Information Table S1) but we have rather considered their geometrical representation showing the spread of models within the 95% confidence interval (Figs. 3a and 3b).

We find (Fig. 3 and Table 2) that all methods based on the unwrapped data better explain the total displacement (ascending+descending) the inversion based on wrapped phases, as shown by the percentage of explained data (87-95% versus 82%) or by the RMS (0.85-2.06×10⁻⁴ versus 2.73×10⁻⁴ m). The dyke inverted from the Conventional inversion is narrower and longer than the initial dyke, resulting in a large distance (94 m) to the initial mesh (Table 2). It is also slightly curved to compensate for the lack of opening resulting from the narrowness of the dyke at depth. These geometrical characteristics are shared by models within the 95% confidence interval (Fig. 3). The narrowness of the dyke induces characteristic residuals on the descending and ascending interferograms with two patterns of positive residuals at the edge of the fracture and a negative patterns to the north-east away from the fracture.

The inversion based on the descending LOS (Conventional Dsc in Fig. 3) leads to a result close to the Conventional inversion conducted with both the ascending and descending interferograms. The dykes obtained are close, being narrow at depth and with an along dip curvature, which are characteristics that were not present in the initial dyke. The average distance between the best model mesh and the initial mesh (61 m) is lower than for the Conventional inversion due to a slightly shorter dyke. However, the use of only one LOS leads to larger uncertainties, a result which is consistent with the study of Dawson & Tregoning (2007) for faults. The residual maps
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Comparison of inverted dyke geometries and residuals for different inversion methods compensating for the data weight (see details in the text). The initial synthetic N140°E-60°E dyke is shown for reference in the top row along with its associated displacement. Panels (a) and (b) give best-fitting model geometries (outlined in black), in map and perpendicular-to-strike views, respectively. Points of colour are mesh nodes from random models generated within the 95% confidence interval of the parameters. (c) shows residuals for descending (left-hand column: Dsc) and ascending (right-hand column: Asc) interferograms with their respective percentage of explained data (%ED) and RMS. See Supporting Information Figs. S3-S7 for comparisons between the subsampled data, the modelled displacements and the residuals.</figcaption></figure>

(Fig. 3c) present the same pattern as the Conventional inversion. This result is in agreement with the descending pass having the most weight in the Conventional inversion. On the other hand, the inversion conducted with displacement projected in an ascending LOS (Conventional Asc) leads to a geometry without any curvature and a fairly constant width, which is closer to the initial dyke. Still, the best-fit model remains longer than the initial dyke, resulting in a high average distance to the initial mesh (97 m). As for the Conventional Descending inversion, the use of the ascending interferogram alone results in large confidence intervals, as indicated by the spread of mesh points within the 95% confidence interval.

The inversion conducted by weighting the covariance matrix leads to a short, straight dyke, which is the closest to the initial mesh (28
<table><caption>**Table 2.** Inversions results for each data weighting method (see inversions method description in text). The percentage of explained data (%ED) and the RMS of the residuals are given for the descending (Dsc), ascending (Asc) and ascending+descending (tot) data, all unmasked points in the interferograms being taken into account. D<sub>Mesh</sub> the average distance between the inverted mesh and the initial mesh and its standard deviation, expressed in meters. (a) gives results for the test on a synthetic N140°E-60°E dyke. (b) concerns the results of the real case of the 2017 July Piton de la Fournaise eruption.</caption><thead><tr><th>Inversion method</th><th>%ED<sub>Dsc</sub></th><th>%ED<sub>Asc</sub></th><th>%ED<sub>tot</sub></th><th>RMS<sub>Dsc</sub> (m)</th><th>RMS<sub>Asc</sub> (m)</th><th>RMS<sub>tot</sub> (m)</th><th>D<sub>Mesh</sub> (m)</th></tr></thead><tbody><tr><td><strong>a. Synthetic test</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Conventional</td><td>95</td><td>80</td><td>92</td><td>1.12×10<sup>-4</sup></td><td>1.26×10<sup>-4</sup></td><td>1.19×10<sup>-4</sup></td><td>94±144</td></tr><tr><td>Conventional Dsc</td><td>94</td><td>78</td><td>91</td><td>1.35×10<sup>-4</sup></td><td>1.39×10<sup>-4</sup></td><td>1.37×10<sup>-4</sup></td><td>61±62</td></tr><tr><td>Conventional Asc</td><td>87</td><td>85</td><td>87</td><td>2.75×10<sup>-4</sup></td><td>0.97×10<sup>-4</sup></td><td>2.06×10<sup>-4</sup></td><td>97±145</td></tr><tr><td>Cov. Weighting</td><td>95</td><td>92</td><td>95</td><td>1.09×10<sup>-4</sup></td><td>0.51×10<sup>-4</sup></td><td>0.85×10<sup>-4</sup></td><td>28±18</td></tr><tr><td> Wrapped</td><td>84</td><td>71</td><td>82</td><td>3.41×10<sup>-4</sup></td><td>1.82×10<sup>-4</sup></td><td>2.73×10<sup>-4</sup></td><td>81±113</td></tr><tr><td><strong>b. 2017 July eruption</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Conventional</td><td>78</td><td>29</td><td>67</td><td>1.37×10<sup>-2</sup></td><td>1.66×10<sup>-2</sup></td><td>1.52×10<sup>-2</sup></td><td>-</td></tr><tr><td>Cov. Weighting</td><td>71</td><td>49</td><td>68</td><td>1.75×10<sup>-2</sup></td><td>1.19×10<sup>-2</sup></td><td>1.50×10<sup>-2</sup></td><td>-</td></tr><tr><td> Wrapped</td><td>78</td><td>23</td><td>66</td><td>1.34×10<sup>-2</sup></td><td>1.79×10<sup>-2</sup></td><td>1.58×10<sup>-2</sup></td><td>-</td></tr></tbody></table>

m). Moreover, the better balance in the ascending and descending passes produces the smallest confidence intervals of all the inversion methods tested here. This method shows the highest percentage of explained data (95%). The goodness of fit is confirmed by the RMS showing the lowest values of all methods for the data set taken individually (1.09×10<sup>-4</sup> m and 0.51×10<sup>-4</sup> m for the descending and ascending, respectively) and in total (0.85×10<sup>-4</sup> m). Such characteristics show that the Covariance Weighting method best accounts for the information contained in both data sets.

The importance of weighting is also indicated by the pressure estimated by each method (Supporting Information Table S1). Indeed, the inversion using only descending data greatly overestimates pressure (0.68 MPa), while the inversion using only ascending data is close (0.53 MPa) to the true value (0.5 MPa). For the Conventional inversion, pressure is intermediate at 0.64 MPa, but, when using the Covariance weighting method, the estimated pressure decreases to 0.57 MPa, which is consistent with a better balance of the data sets weights.

Finally, the Wrapped inversion, shows a dyke geometry similar to the Conventional inversion. This dyke tends to deviate eastwards at depth. The confidence intervals associated with this inversion are slightly lower than those for the Conventional inversion. The percentage of explained data is the lowest (82%), and the residuals map shows high magnitude residuals which can be attributed to the poor fit of the data magnitude.

## 4 Application to the 2017 July Eruption of Piton de la Fournaise

Piton de la Fournaise is a hotspot volcano located on La Réunion Island (Indian Ocean) (Fig. 4 for location). Eruptions at Piton de la Fournaise are typically fissure eruptions. In more than 70% of the cases, they occur along two rift zones located to the north-east (NE) and the south-east (SE) of the summit Dolomieu crater (Fig. 4 and Michon et al. (2007)). They may also occur along a N120° rift zone corresponding to an axis linking the extinct Piton des Neiges volcano and Piton de la Fournaise.

The co-eruptive displacements recorded by InSAR show that dyke intrusions produce displacement with a clear asymmetrical pattern characterized by large displacements east of the eruptive fissures compared to the west (Peltier 2007; Fukushima et al. 2010; Richter & Froger 2020). This characteristic is shown both on the cumulated displacements from October 2010 to April 2020 (Supporting Information Fig. S8), as well as on the displacement patterns of individual rift zone eruptions, whether they took place in the north, middle or south rift zones (Supporting Information Fig. S9). Such an asymmetrical displacement produces significantly different patterns for InSAR data acquired in ascending and descending orbital passes. Typically, descending data commonly show one main displacement pattern to the east of the eruptive fissure, corresponding to a ground motion towards the satellite, whereas ascending data show a displacement pattern exhibiting two lobes with one corresponding to a ground motion towards the satellite and one away from it. Both lobes have a lower magnitude than the pattern of the descending pass (Fig. 5 and Supporting Information Fig. S8).

The 2017 July eruption was a typical SE Rift Zone eruption. The displacement induced by the dyke intrusion was imaged on four independent interferograms: two pairs of ascending and descending interferograms from Sentinel-1 and COSMO-Skymed satellites (Supporting Information Table S2). We inverted this eruption using the three inversion approaches presented above with the same search intervals (Supporting Information Table S3). We only inverted the ascending/descending pair from Sentinel-1 because of their better quality, although COSMO-Skymed images are used to discuss possible sources of residuals.

Before being subsampled, interferograms were unwrapped using the SNAPHU algorithm (Chen & Zebker 2000, 2001) within an iterative procedure developed by Fukushima et al. (2005). As Piton de la Fournaise area is small (maximum of 20x20 km), meteorological models
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Shaded relief map of Piton de la Fournaise volcano with main locations. (a) Island location in Indian ocean. (b) Réunion Island. (c) Location area of the study. Green dashed square indicates the area shown in Figs. 5-7 and Supporting Information Figs. S8, S9, S12, S13 and S14. Red lines show eruptive fissure locations since 1998. RZ: rift zones.</figcaption></figure>

are too coarse and are therefore not relevant to correct atmospheric effects. Consequently, we correct "large scale" (km to tens km) atmospheric effects by filtering out large wavelength signal through the subtraction of the best-fitting polynomial surface. This method works well at Piton de la Fournaise as eruption-induced deformation only affect part of the Enclos Fouqué caldera, so that polynomial surface can be adjusted using the rest of the Enclos Fouqué. As our study area is very small, LOS values vary by $6 \times 10^{-3}$ within each InSAR image despite the elevation change which is negligible with respect to measurements and modelling standard deviations (4 µm of differences for a displacement of 1 cm in EW, NS and UD). We therefore use a constant LOS value for ascending and descending data, respectively.

The sub-sampling was carried out independently on ascending and descending interferograms after masking points affected by the lava flow (magenta contour on Fig. 5a), as well as points within Dolomieu crater and highly incoherent pixels. Note that points outside the caldera rim are generally highly incoherent due to the dense vegetation and are therefore masked. Similarly as for the synthetic tests, we have adjusted the parameters of the quadtree algorithm in order to obtain almost the same number of points for the ascending and descending data sets (780 and 776 points, respectively). Both Sentinel-1 interferograms and the subsampled data show significant differences in magnitudes (Fig. 5). Descending data reach values higher than 0.4 m, whereas the ascending data are lower than 0.1 m leading to a relative weight of the descending over total displacement $W_{Dsc} \approx 0.86$.

Similarly as in section 2.1, the initial covariance matrix was computed with a noise variance $\sigma^2 = 5 \times 10^{-4} \text{ m}^2$ and a correlation distance of 850 m estimated by Fukushima et al. (2010) for the Piton de la Fournaise. In the case of the Covariance Weighting inversion, we took the initial covariance matrix and we weighted it using eq. (5). For the Wrapped inversion, the decrease of spatial resolution caused by subsampling could result in a loss of the phase continuity (aliasing). Therefore, we took the unwrapped interferograms and then rewrapped data with a wavelength large enough to avoid any aliasing in the sub-sampled data. By trials and errors, we choose an half-wavelength of 8 cm to make sure the sampling points were appropriately distributed in each fringe. The covariance matrix used the same variance as before ($5 \times 10^{-4} \text{ m}^2$), but converted to phase values.

The total percentage of explained data is very close for all three methods (66-68%) as well as the final total RMS between 1.50 and 1.58 cm (Table 2). Variance of the residuals goes up to $7 \times 10^{-4} \text{ m}^2$ which is very close to the value $(5 \times 10^{-4} \text{ m}^2)$ considered for this study. The best-fitting model of the Conventional inversion is able to explain 67% of the data. In details, the descending interferograms is explained at 78% but the ascending one is only explained at 29%. Similarly, the Wrapped inversion better explains the descending (78%)
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Ascending and descending sentinel-1 data for the 2017 July eruption of Piton de la Fournaise. (a) Unwrapped interferograms in Descending (left-hand panels) and Ascending (right-hand panels) passes. Top panel gives the full data while the bottom panel shows the interferograms subsampled by a quadtree algorithm. (b) Subsampled displacement values given as a data vector in an arbitrary order. Note the clear magnitude difference between the descending (Dsc) and ascending (Asc) data sets, with a ratio of ≈ 5, which leads to a relative weight $W_{Dsc}$ ≈0.86.</figcaption></figure>

than the ascending (23%) interferogram. For both methods, RMS shows lower values for the descending data set despite its higher average magnitude, and higher values for the ascending data set, indicating that data sets fit is independent of their magnitude. In contrast, the Covariance Weighting inversion is 20% better at explaining the ascending interferogram than the Conventional inversion, at the expanse of a decrease of 7% in the percentage of explained data for the descending interferogram. Similarly, the RMS increases for the descending data set and decreases for the ascending one, leading to the descending data set having a higher RMS than the ascending one, consistently with their relative average magnitudes. In terms of inverted dyke geometries (Fig. 6), the three inversions show similar characteristics: an elongated flat part which then becomes vertical close to the eruptive fissure. The main differences between the three inversion methods are: (1) the lack of bottom line curvatures for the Covariance weighting and the Wrapped inversions; (2) the dyke maximum opening being more centered beneath the summit using the Covariance Weighting inversion; and (3) the Weighted and Wrapped inversions show clearer maximum in the marginal posterior probability density function than the Conventional inversion (Supporting Information Fig. S11), indicating that the best-fit model is better determined. This reflects in the spread of mesh points corresponding to the 95% confidence interval (Fig. 7).

Residual maps (Fig. 6) show that the elongated dyke determined by the Weighted inversion best explains the ascending data. With the Conventional inversion, the descending displacements are overestimated, while the ascending ones are underestimated. On the other hand, inversion methods designed to give a better balance between data, underestimate displacements on the east side of the Dolomieu crater and overestimate displacements on the west side of the crater both on the descending and the ascending passes. This indicates that a better compromise is achieved. Globally, the percentage of explained data is similar for all the methods (Table 2), but the relative percentage of explained data for each data sets is more balanced using the Covariance Weighting inversion.

# 5 Discussion

The Wrapped inversion, successfully fits both the ascending and descending interferograms. It produces results similar to the Conventional inversion for the synthetic tests and the Piton de la Fournaise July 2017 eruption, confirming the relevance of its use when conducting inversions without any unwrapping. However the signal magnitude is poorly fitted in the synthetic tests, due to overly low pressure, while it is well fitted in the inversion of the July 2017 eruption. We attribute this inexact pressure determination to the subsampling strategy. In order, to be able to constrain the source pressure, the displacement gradients must be accurately accounted for
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Best-fitting model obtained for the July 2017 eruption of Piton de la Fournaise using three different inversion methods: Conventional, Covariance Weighting and Wrapped (see text for detailed description). (a) Wrapped Sentinel-1 descending (Dsc) and ascending (Asc) interferograms. (b) Modelled and residual displacements for each method. Green line represents the eruptive fissure and pink outline represents the map view of the intrusion contour. (c) Corresponding 3-D best-fitting model.</figcaption></figure>

and all fringes must be represented in the subsampled data. The specific sub-sampling strategy for wrapped data as presented recently by Jiang & González (2020) could address this issue.

In the example of the 2017 Piton de la Fournaise eruption, giving more weight to the ascending LOS by using the Covariance Weighting inversion results in an intrusion rooted beneath Dolomieu crater. This geometry is in better agreement with pre-eruptive seismicity, which is located beneath Dolomieu (Lengliné et al. 2016). Previous studies from InSAR data on Piton de la Fournaise noted that intrusions have to root to the east of Dolomieu crater (Froger et al. 2004; Fukushima et al. 2005, 2010) to account for the strong displacement asymmetry between the west and east sides of the eruptive fissures (Supporting Information Figs. S8 and S9). As an alternative to these models, Got et al. (2013) proposed that a plastic behaviour in the shallowest part of the volcano edifice could also explain the displacement asymmetry. However, their model was only two-dimensional and only considers GNSS data to compare model and observations. Consequently, there is no evidence that their model could be acceptable for 3-D displacements. In our study, the use of curved dykes together with a better balance between the ascending and descending data sets leads to intrusions that root beneath Dolomieu crater, rather than to the east, with no need for inelastic rock behaviour.

It is always possible to obtain a better fit by adding complexity to a model. However, a better fit does not necessarily mean that a model
<figure><img src="image_8.png" /><figcaption>Figure 7. Inversion of the July 2017 eruption. Mesh points corresponding to the 95% confidence interval for the Conventional, Covariance Weighting and Wrapped inversions are represented. Colour-scale represents the normalized probability density of mesh nodes.</figcaption></figure>

is better. In order to test that adding curvatures leads to better models, we rely on the Akaike Information Criteria (Akaike 1974):

$$AIC = 2k + \chi^2 + cst \quad (9)$$

where k is the number of model parameters, $\chi^2$ is the misfit and cst is a constant that depends on the number of subsampling points and the covariance matrix. The better a model, the lower the criteria. We found that the AIC is lower when considering the curvature parameters, indicating that models with curvatures are preferable (see Supporting Information Fig. S12). Moreover, thanks to the curvature parameters, the near surface dip of the intrusion is sub-vertical which is consistent with field data showing mostly vertical dykes (Le-tourneur et al. 2008). Whereas, without curvatures, the intrusion has a much lower dip (≈50°). The curvature leads to a transition from a sill at depth to a dyke shallower in the edifice. Such transition was already found at Piton de la Fournaise for the May 2016 eruption (Smittarello et al. 2019b). It was also observed at Fernandina Volcano in Galápagos in the case of circumferential intrusion (Chadwick et al. 2011; Bagnardi et al. 2013). Theoretical occurrence of sill to dyke transitions was shown by numerical modelling of Fialko (2001) for intrusion in mid-ocean ridge context. Therefore, we consider this sill to dyke geometry to be likely.

After correcting for unbalanced magnitudes, data sets from ascending and descending orbits still have 20% of difference in their percentage of explained data (Table 2), with displacements to the east of the intrusion being underestimated and displacements to the west of the intrusion being overestimated. Residuals could originate from (1) a supplementary source, (2) interferograms covering different time periods or (3) an incorrect or too simplistic forward model. Here we discuss these hypothesis. (1) The presence of supplementary sources is inconsistent with the data as discussed by Smittarello et al. (2019a). (2) InSAR and GNSS data show no noticeable deformation before the ascending acquisition or after the descending acquisition, indicating that the observed residuals can not come from variable time spans of the interferometric pairs (Fig. 8, and interferograms given in Supporting Information Fig. S13). Despite an acquisition delay of 8 days (Fig. 8), the COSMO-Skymed interferograms record the same deformation, and show the same residual pattern as our model (Supporting Information Fig. S14), indicating that the ascending and descending data sets recorded the same deformation. This is in agreement with the typical behaviour of Piton de la Fournaise where intrusive events are accompanied by rapid deformation without any significant pre-eruptive or post-eruptive deformation (Peltier 2007; Fukushima et al. 2010). (3) Finally, model assumptions, such as the source type, source description, or the mechanical properties of the host medium could be too simple. The intrusion geometry might be more complex than a quadrangle and the driving pressure might not be constant. In situ observations of dykes show that intrusions do not necessarily follow a planar surface (Delaney et al. 1986; Maccaferri et al. 2011; Kavanagh et al. 2006; Menand 2008) and that openings do not simply results from a constant overpressure (Kavanagh & Sparks 2011).

## 6 Conclusion

By means of synthetic tests, we show that the LOS of ascending and descending data sets influences the weight of each subsampled data set considered in an inversion. For roughly 60° dipping dykes striking NS, the data set with the largest magnitude accounts for 90% of the data, resulting in a biased geometry and pressure of the inverted model.

To mitigate the influence of the relative magnitudes of a data set, we compared a Conventional inversion with two methods: a Covariance Weighting approach where the data covariance matrix is tuned to give the same weight to each data set, and a Wrapped inversion approach where the inversion is conducted on the wrapped data. Synthetic tests show that the Covariance Weighting inversion leads to a more balanced fit between data sets than a Conventional inversion. More importantly, these tests show that the final model geometry is
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Available interferograms for the July 2017 eruption. Cyan and blue bars show time spans of ascending and descending interferograms, respectively, for COSMO-Skymed (CSK) and Sentinel-1 (S1) satellites. The red band gives eruption period and the grey bar represents the possible time window for deformation to have occurred.</figcaption></figure>

closer to the initial model and has smaller confidence intervals. In contrast, the Wrapped Data inversion leads neither to a more balanced data fit nor to a better model than the Conventional inversion.

For the July 2017 eruption of Piton de la Fournaise, the Covariance Weighting inversion improves the data fit balance. The difference in the percentage of explained data for the descending and ascending data sets is reduced from 50% to 20%. This method also leads to a slightly different geometry, with an intrusion rooted beneath the summit crater where pre-eruptive seismicity is usually recorded, rather than to the east of this crater as for the Conventional inversion and previous studies at this volcano.

Our study concerns dykes, but the consequence of unbalanced magnitude between ascending and descending InSAR data sets may also affect other sources of asymmetric displacements, such as faults, for which we recommend using our Covariance Weighting method. Moreover, the weighting method developed for GNSS and InSAR data (Yabuki & Matsu'ura 1992; Tarantola 2005; Fukuda & Johnson 2010; Auriac et al. 2013; Sigmundsson et al. 2014) can be used jointly with our method in order to control the relative weighting between GNSS and InSAR data, but also between ascending and descending InSAR data sets.

## ACKNOWLEDGMENTS

This research was financed by the French Government Laboratory of Excellence initiative n° ANR-10-LABX-0006, the Region Auvergne and the European Regional Development Fund. This is Laboratory of Excellence ClerVolc contribution number 473. All interferograms were produced with the Diapason software and inversions have been performed on the supercomputer facilities of the Mésocentre Clermont-Auvergne. We are thankful to Delphine Smittarello for constructive discussion on data weighting and for proofreading this manuscript. Finally, we thank Frances van Wyk de Vries for proofreading the English of this manuscript.

## DATA AVAILABILITY

InSAR data were acquired within the framework of the Indian Ocean InSAR Observatory (OI²/OPGC/SNOV/INSU), where the interferograms used in this study are available on the CASOAR link (https://wwwobs.univ-bpclermont.fr/casoar). Otherwise, unwrapped interferograms, meshes and results of simulations can be downloaded via the Zenodo website (https://doi.org/10.5281/zenodo.4629171). To launch inversions presented in this study, the DEFVOLC software is available online on the OPGC website (http://opgc.fr/defvolc).

## REFERENCES

* Akaike, H., 1974. A New Look at the Statistical Model Identification, *IEEE Transactions on Automatic Control*, **19**(6), 716-723.

* Ali, S. T. & Feigl, K. L., 2012. A new strategy for estimating geophysical parameters from InSAR data: Application to the Krafla central volcano in Iceland, *Geochemistry, Geophysics, Geosystems*, **13**(6), 1-18.

* Árnadóttir, T., Geirsson, H., & Einarsson, P., 2004. Coseismic stress changes and crustal deformation on the Reykjanes Peninsula due to triggered earthquakes on 17 June 2000, *Journal of Geophysical Research: Solid Earth*, **109**(9), 1-12.

* Auriac, A., Spaans, K. H., Sigmundsson, F., Hooper, A., Schmidt, P., & Lund, B., 2013. Iceland rising: Solid Earth response to ice retreat inferred from satellite radar interferometry and visocelastic modeling, *Journal of Geophysical Research: Solid Earth*, **118**(4), 1331-1344.

* Bagnardi, M., Amelung, F., & Poland, M. P., 2013. A new model for the growth of basaltic shields based on deformation of Fernandina volcano, Galápagos Islands, *Earth and Planetary Science Letters*, **377-378**, 358-366.

* Carrier, A., Got, J.-l., Peltier, A., Ferrazzini, V., Staudacher, T., Kowalski, P., & Boissier, P., 2015. A damage model for volcanic edifices: Implications for edifice strength, magma pressure, and eruptive processes, *Journal of Geophysical Research: Solid Earth*, pp. 1-17.

* Cayol, V. & Cornet, F. H., 1997. 3D mixed boundary elements for elastostatic deformation field analysis, *International journal of rock mechanics and mining*..., **34**(2), 275-287.
Cayol, V. & Cornet, F. H., 1998. Three-dimensional modeling of the 1983–1984 eruption at Piton de la Fournaise Volcano, Réunion Island, *Journal of Geophysical Research*, **103**(B8), 18025.

Cayol, V., Fukushima, Y., Dabaghi, F., Bodart, O., Tridon, M., Smitarello, D., Dumont, Q., Froger, J.-L., Wauthier, C., & Cattelain, Cyprien Gouinaux, C., 2020. DefVolc: tools for 3D inverse modelling of deformation data.

Chadwick, W. W., Jónsson, S., Geist, D. J., Poland, M., Johnson, D. J., Batt, S., Harpp, K. S., & Ruiz, A., 2011. The May 2005 eruption of Fernandina volcano, Galápagos: The first circumferential dike intrusion observed by GPS and InSAR, *Bulletin of Volcanology*, **73**(6), 679–697.

Chen, C. W. & Zebker, H. A., 2000. Two-Dimensional Phase Unwrapping with Statistical Models for Nonlinear Optimization, *IEEE*.

Chen, C. W. & Zebker, H. A., 2001. Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization, *Journal of the Optical Society of America A*, **18**(2), 338.

Chen, Y., Remy, D., Froger, J. L., Peltier, A., Villeneuve, N., Darrozes, J., Perfettini, H., & Bonvalot, S., 2017. Long-term ground displacement observations using InSAR and GNSS at Piton de la Fournaise volcano between 2009 and 2014, *Remote Sensing of Environment*, **194**(March 1998), 230–247.

Dawson, J. & Tregoning, P., 2007. Uncertainty analysis of earthquake source parameters determined from InSAR: A simulation study, *Journal of Geophysical Research: Solid Earth*, **112**(9), 1–13.

Delaney, P. T., Pollard, D. D., Ziony, I., & Mckee, E. H., 1986. Field Relations Between Dikes and Joints: Emplacement Processes and Paleostress Analysis, *JOURNAL OF GEOPHYSICAL RESEARCH*, **91**(5), 4920–4938.

Doin, M. P., Lasserre, C., Peltzer, G., Cavalié, O., & Doubre, C., 2009. Corrections of stratified tropospheric delays in SAR interferometry: Validation with global atmospheric models, *Journal of Applied Geophysics*, **69**(1), 35–50.

Duputel, Z., Jiang, J., Jolivet, R., Simons, M., Rivera, L., Ampuero, J. P., Riel, B., Owen, S. E., Moore, A. W., Samsonov, S. V., Ortega Culaciati, F., & Minson, S. E., 2015. The Iquique earthquake sequence of April 2014: Bayesian modeling accounting for prediction uncertainty, *Geophysical Research Letters*, **42**(19), 7949–7957.

Feigl, K. L. & Thurber, C. H., 2009. A method for modelling radar interferograms without phase unwrapping: Application to the M 5 Fawnskin, California earthquake of 1992 December 4, *Geophysical Journal International*, **176**(2), 491–504.

Fialko, Y., 2001. On origin of near-axis volcanism and faulting at fast spreading mid-ocean ridges, *Earth and Planetary Science Letters*, **190**(1-2), 31–39.

Fialko, Y., 2004. Probing the mechanical properties of seismically active crust with space geodesy: Study of the coseismic deformation due to the 1992 Mw 7.3 Landers (southern California) earthquake, *Journal of Geophysical Research*, **109**(B3), 1–19.

Fialko, Y., Simons, M., & Agnew, D., 2001. The complete (3-D) surface displacement field in the epicentral area of the 1999 M7.1 Hector Mine earthquake California from space geodetic observations, *Geophys. Res. Lett.*, **28**(16), 3063–3066.

Froger, J. L., Fukushima, Y., Briole, P., Staudacher, T., Souriot, T., & Villeneuve, N., 2004. The deformation field of the August 2003 eruption at Piton de la Fournaise, Reunion Island, mapped by ASAR interferometry, *Geophysical Research Letters*, **31**(14), 1–5.

Fukuda, J. & Johnson, K. M., 2008. A fully Bayesian inversion for spatial distribution of fault slip with objective smoothing, *Bulletin of the Seismological Society of America*, **98**(3), 1128–1146.

Fukuda, J. & Johnson, K. M., 2010. Mixed linear-non-linear inversion of crustal deformation data: Bayesian inference of model, weighting and regularization parameters, *Geophysical Journal International*, **181**(3), 1441–1458.

Fukushima, Y., Cayol, V., & Durand, P., 2005. Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise, *Journal of Geophysical Research B: Solid Earth*, **110**(3), 1–15.

Fukushima, Y., Cayol, V., Durand, P., & Massonnet, D., 2010. Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island, *Journal of Geophysical Research: Solid Earth*, **115**(10).

Got, J. L., Peltier, A., Staudacher, T., Kowalski, P., & Boissier, P., 2013. Edifice strength and magma transfer modulation at Piton de la Fournaise volcano, *Journal of Geophysical Research: Solid Earth*, **118**(9), 5040–5057.

Heap, M. J., Villeneuve, M., Albino, F., Farquharson, J. I., Brothelande, E., Amelung, F., Got, J. L., & Baud, P., 2020. Towards more realistic values of elastic moduli for volcano modelling, *Journal of Volcanology and Geothermal Research*, **390**, 106684.

Jiang, Y. & González, P. J., 2020. Bayesian Inversion of Wrapped Satellite Interferometric Phase to Estimate Fault and Volcano Surface Ground Deformation Models, *Journal of Geophysical Research: Solid Earth*, **125**(5).

Jónsson, S., Zebker, H., Segall, P., & Amelung, F., 2002. Fault Slip Distribution of the 1999 Mw 7.1 Hector Mine, California, Earthquake, Estimated from Satellite Radar and GPS Measurements, *Bulletin of the Seismological Society of America*, **92**(41), 1377–1389.

Kavanagh, J. L. & Sparks, R. S. J., 2011. Insights of dyke emplacement mechanics from detailed 3D dyke thickness datasets, *Journal of the Geological Society*, **168**, 965–978.

Kavanagh, J. L., Menand, T., & Sparks, R. S. J., 2006. An experimental investigation of sill formation and propagation in layered elastic media, *Earth and Planetary Science Letters*, **245**(3-4), 799–813.

Lengliné, O., Duputel, Z., & Ferrazzini, V., 2016. Uncovering the hidden signature of a magmatic recharge at Piton de la Fournaise volcano using small earthquakes, *Geophysical Research Letters*, **43**(9), 4255–4262.

Letourneur, L., Peltier, A., Staudacher, T., & Gudmundsson, A., 2008. The effects of rock heterogeneities on dyke paths and asymmetric ground deformation: The example of Piton de la Fournaise (Réunion Island), *Journal of Geophysical Research: Solid Earth*, **112**(4), 1–21.

Lohman, R. B. & Simons, M., 2005. Some thoughts on the use of InSAR data to constrain models of surface deformation: Noise structure and data downsampling, *Geochemistry, Geophysics, Geosystems*, **6**(1).

Maccaferri, F., Bonafede, M., & Rivalta, E., 2011. A quantitative study of the mechanisms governing dike propagation, dike arrest and sill formation, *Journal of Volcanology and Geothermal Research*, **208**(1-2), 39–50.

Menand, T., 2008. The mechanics and dynamics of sills in layered elastic rocks and their implications for the growth of laccoliths and other igneous complexes, *Earth and Planetary Science Letters*, **267**(1-2), 93–99.

Michon, L., Saint-Ange, F., Bachelery, P., Villeneuve, N., & Staudacher, T., 2007. Role of the structural inheritance of the oceanic lithosphere in the magmato-tectonic evolution of Piton de la Fournaise volcano (La Réunion Island), *Journal of Geophysical Research: Solid Earth*, **112**(4), 1–21.

Montgomery-Brown, E. K., Sinnett, D. K., Larson, K. M., Poland, M. P., Segall, P., & Miklius, A., 2011. Spatiotemporal evolution of dike opening and décollement slip at Kilauea Volcano, Hawai'i, *Journal of Geophysical Research: Solid Earth*, **116**(3), 1–14.

Pascal, K., Neuberg, J., & Rivalta, E., 2013. On precisely modelling surface deformation due to interacting magma chambers and dykes, *Geophysical Journal International*, **196**(1), 253–278.

Peltier, A., 2007. Suivi, modélisation et évolution des processus d'injections magmatiques au Piton de La Fournaise (Réunion) à partir d'une analyse croisée des données de déformation, géochimiques et structurales, *Géologie appliquée. Université de la Réunion*.

Peltier, A., Froger, J. L., Villeneuve, N., & Catry, T., 2017. Assessing the reliability and consistency of InSAR and GNSS data for retrieving 3D-displacement
rapid changes, the example of the 2015 Piton de la Fournaise eruptions, *Journal of Volcanology and Geothermal Research*, **344**, 106–120.

Remy, D., Froger, J.-L., Perfettini, H., Bonvalot, S., Gabalda, G., Albino, F., Cayol, V., Legrand, D., & De Saint Blanquat, M., 2014. Persistent uplift of the Lazufre volcanic complex (Central Andes): New insights from PCAIM inversion of InSAR time series and GPS data, *Geochemistry, Geophysics, Geosystems*, **15**, 4692–4711.

Richter, N. & Froger, J. L., 2020. The role of Interferometric synthetic aperture radar in detecting, mapping, monitoring, and modelling the volcanic activity of piton de la Fournaise, La Reunion: A review, *Remote Sensing*, **12**(6).

Sambridge, M., 1998. Exploring multidimensional surfaces without a map, *Inverse Problems*, **14**(3), 427–440.

Sambridge, M., 1999a. Geophysical inversion with a neighbourhood algorithm - 1. Searching a parameter space, *Geophysical Journal International*, **138**, 479–494.

Sambridge, M., 1999b. Geophysical inversion with a neighbourhood algorithm - 2. Appraising the ensemble, *Geophysical Journal International*, **138**, 727–746.

Sigmundsson, F., Hooper, A., Hreinsdóttir, S., Vogfjörd, K. S., Ófeigsson, B. G., Heimisson, E. R., Dumont, S., Parks, M., Spaans, K., Gudmundsson, G. B., Drouin, V., Árnadóttir, T., Jónsdóttir, K., Gudmundsson, M. T., Högnadóttir, T., Fridriksdóttir, H. M., Hensch, M., Einarsson, P., Magnússon, E., Samsonov, S., Brandsdóttir, B., White, R. S., Ágústsdóttir, T., Greenfield, T., Green, R. G., Hjartardóttir, Á. R., Pedersen, R., Bennett, R. A., Geirsson, H., la Femina, P. C., Björnsson, H., Pálsson, F., Sturkell, E., Bean, C. J., Möllhoff, M., Braiden, A. K., & Eibl, E. P., 2014. Segmented lateral dyke growth in a rifting event at Bárðarbunga volcanic system, Iceland, *Nature*, **517**(7533).

Smittarello, D., Cayol, V., Pinel, V., Froger, J. L., Peltier, A., & Dumont, Q., 2019a. Combining InSAR and GNSS to track magma transport at basaltic volcanoes, *Remote Sensing*, **11**(19), 1–26.

Smittarello, D., Cayol, V., Pinel, V., Peltier, A., Froger, J.-L., & Ferrazzini, V., 2019b. Magma Propagation at Piton de la Fournaise From Joint Inversion of InSAR and GNSS, *Journal of Geophysical Research: Solid Earth*, **124**(2), 1361–1387.

Sudhaus, H. & Jónsson, S., 2009. Improved source modelling through combined use of InSAR and GPS under consideration of correlated data errors: Application to the June 2000 Kleifarvatn earthquake, Iceland, *Geophysical Journal International*, **176**(2), 389–404.

Tarantola, A., 1987. *Inverse Problem Theory*, Elsevier.

Tarantola, A., 2005. *Inverse Problem Theory*, vol. 120, Society for Industrial and Applied Mathematics.

Tridon, M., Cayol, V., Froger, J. L., Augier, A., & Bachèlery, P., 2016. Inversion of coeval shear and normal stress of Piton de la Fournaise flank displacement, *Journal of Geophysical Research: Solid Earth*, **121**(11), 7846–7866.

Welstead, S. T., 1999. *Fractal and wavelet image compression techniques*, no. TT 40, The International Society for Optical Engineering.

Wright, T. J., Parsons, B. E., & Zhong, L., 2004. Toward mapping surface deformation in three dimensions using InSAR, *Geophysical Research Letters*, **31**(1), L01607.

Yabuki, T. & Matsu'ura, M., 1992. Geodetic data inversion using a Bayesian information criterion for spatial distribution of fault slip, *Geophysical Journal International*, **109**(2), 363–375.

Zebker, H. A., Rosen, P. A., & Hensley, S., 1997. Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps, *Journal of Geophysical Research*, **102**, 7547–7563.