*3.1. SO2 Column Density Retrieval*

SO2 path-length concentrations were measured in the plume of Láscar volcano using a NOVAC-type stationary scanning Mini-DOAS instrument (Figure 1; [49]; NOVAC: Network for the Observation of Volcanic and Atmospheric Change) that was permanently installed in April 2013 by the Chilean Observatorio Volcanológico De los Andes del Sur (OVDAS). The DOAS station (Lejía DOAS) was deployed east of Láscar (at Lat −23.387◦, Long −67.678◦), according to prevailing westerly wind directions. The DOAS-instrument scans across the sky from horizon to horizon along a semi-conical surface (Figure 2a) and measures the spectra of the incoming scattered sunlight at 51 angular steps of 3.6◦ [49]. Each complete scan takes about 5–10 min and yields a crosswind SO2 total column profile, i.e., perpendicular to the transport direction of the volcanic gas plume. Measurements are conducted during daylight, i.e., when UV intensities are sufficiently high to achieve an acceptable signal-to-noise ratio. At these latitudes the required conditions are met during a time interval approximately ranging from about 11:00 a.m. to 10:00 p.m. (UTC) in austral summer, and 12:00 a.m. to 09:00 p.m. (UTC) in austral winter.

SO2 differential slant column densities (SCDs) were retrieved from the sunlight spectra by means of the DOAS method [50], as implemented in an automated evaluation routine of the NOVAC-software [49]. Evaluation was performed in the wavelength range 310–325 nm using absorption cross sections of SO2 [51], and O3 [52], which were convolved to the spectral resolution of the instrument. In addition, a Ring spectrum was included in the DOAS-fit, in order to avoid the filling in of the Fraunhofer lines of the solar spectrum as a result of rotational Raman scattering [53]. Shift correction was applied to the measured spectra using the known features of a pixel-wavelength calibrated solar spectrum, which was derived from the Solar flux atlas of Kurucz et al. [54].

**Figure 2.** Spatio-temporal relationships between the methods used for the gas plume estimate. (**a**) Sketch showing the measurement geometry and location of the Multi-GAS and scanning Mini-DOAS instruments at Láscar volcano during an overpass of TerraSAR-X. Refractive delay of the radar occurs inside the volcanic plume, which is heading towards East, due to predominantly westerly winds. See text for discussion; (**b**) Date versus time plot depicting acquisition times of scanning DOAS measurements and SAR images. Measurement times are indicated using Coordinated Universal Time (UTC), which is offset by +3 h with respect to Chile Summer Time (CLST). *Red arrows* indicate the temporal offset between SAR images and scanning DOAS data chosen for analysis; (**c**) Spatial and temporal baselines of SAR images used in DInSAR time series subsets 01 (*blue*) and 02 (*red*). Master scenes of *subset 01* and *subset 02* are from 18 October 2013 and 12 December 2013, respectively, and each computed interferogram is represented by a line between the corresponding two images.

#### *3.2. Estimation of Water Vapor Contents in the Láscar Plume*

Water vapor contents in the Láscar gas plume were obtained by scaling the SO2 SCDs with the molar H2O/SO2 ratio determined from gas concentration measurements, which were in turn conducted using a multi-component Gas Analyzer System [55] during a field survey on 2 December 2012 (see [33] for details), that is roughly one year in advance of the period that we consider here for the analysis of our SAR and Mini-DOAS observations (Figure 2b). The Multi-GAS instrument basically consists of an air pump feeding ambient air into a suite of electrochemical sensors for the measurement of SO2, H2S and H2, and a non-dispersive closed-path IR spectrometer for measuring CO2 and H2O concentrations. Gas concentrations were measured for several hours in a dilute and partly condensed volcanic plume crossing the southern rim of the active crater (roughly at Lat −23.366◦, Long −67.734◦; see Figure 2a), and the signals of all sensors were recorded at 0.5 Hz by a data-logger, yielding a large number of measurements. The raw data was processed according to [56,57], in order to identify the characteristic molar gas ratios [33].

Multi-GAS instruments, however, exclusively measure gaseous water vapor and do not take into account the liquid water content (LWC) in the volcanic cloud. Due to the obvious condensed nature of the plume at the crater rim, our PWV contents were thus estimated on the basis of the maximum molar H2O/SO2 ratio obtained for the measurement period, which was used to scale the SO2 SCDs yielding a first PWV estimate. This first estimate, we however consider to be representative solely for the proximal portion of the volcanic plume (see Appendix A for details on the conversion of SO2 SCDs to PWV contents).

To estimate PWV contents also from the distal portion of the volcanic plume, we additionally consider condensation and evaporation processes, which change the liquid-to-vapor ratio during transport of the volcanic plume. LWCs are typically elevated close to the emission source and decrease with increasing downwind distance, due to gradual entrainment of dry ambient air, causing evaporation of cloud droplets (e.g., [58]). This in turn results in an increase of water vapor in the downwind portion of the volcanic cloud, which may further be promoted by the presence of aerosols [59], and which is expected to be reflected as an enhancement of the radar delay in interferometric measurements.

Potential evaporation rates were thus determined for the measurement periods of the scanning DOAS and then used to upscale our fixed molar gas ratio, which was thereby adjusted to increase proportionally to evaporation rates, which on the other hand strongly depend on wind speeds over arid areas (see Appendix B for details on evaporation calculation). The resulting variable ratio was used to upscale the SO2 SCDs at each scan angle of the plume cross sections as described above (and in Appendix A), yielding a second PWV estimate, which we deem to be representative for the distal part of the gas plume.

#### *3.3. SAR Data and InSAR Methods*

The aforementioned water vapor maps into the used DInSAR data, which are now briefly described. TerraSAR-X was tasked to acquire high resolution spotlight SAR scenes covering a 10 × 10 km large area around Láscar volcano on both ascending and descending passes of the orbit. TerraSAR-X flies along the day-night boundary of the Earth in a near-polar sun-synchronous dusk-dawn orbit with an 11-day repeat cycle. SAR observations of a given point on Earth's surface are thus conducted within fixed time windows, which at Láscar correspond to day times of about 10:00 a.m. (UTC) on descending nodes and 11:00 p.m. (UTC) on ascending nodes of the orbit. The DInSAR time series considered here covers a 5 months period from 18 October 2013 to 16 February 2014 (Figure 2b). DInSAR preparation and analysis was done using the dedicated modular InSAR software system GENESIS from DLR [60]. A set of 9 SAR images from descending track 111 was chosen according to availability of complementary scanning DOAS data (Figure 2b), and these images were used to create two temporally overlapping disjoint interferogram time series, termed *subset 01* and *subset 02* hereafter. The subsets consist of 3 and 4 temporally interconnected interferograms, respectively, i.e.,

the interferograms of each subset are based on a common master scene (Figure 2c; Table 1). SAR images were combined to form interferograms with very small spatial baselines, which accordingly resulted in different temporal baselines for each interferometric pair as depicted in Figure 2c and Table 1. The very small spatial baselines chosen here minimize phase contributions from topography and prevent from unwrapping failures due to a coarse digital elevation model (DEM).



#### *3.4. Determination of Gas Plume Related Phase Delays*

PWV contents determined for the gas plumes of each SAR acquisition were used to compute the corresponding differential phase delays encountered in each interferogram. For this purpose we first determined the slant wet delay (SWD) that such PWV contents theoretically would produce in the LOS of a SAR image. Determination of the theoretical SWD is straightforward, since it can be inferred from the zenith path delay (ZPD), which is linearly related to the vertically integrated water vapor (IWV) contents, where 1 mm in PWV contents roughly results in a 6.5 mm ZPD of microwave signals [61,62], and the ZPD accordingly is

$$
\delta^{ZP} = II\_{T\_S}^{-1} \text{ I\!\!\!W\!V\!V\!/} \tag{1}
$$

where *<sup>δ</sup>ZP* is the ZPD in (mm), *IWV* is the columnar integrated water vapor in (kg·m−2), which is approximately equal to PWV contents in (mm), and *Π*−<sup>1</sup> *TS* (≈ 6.5) is a dimensionless conversion factor that is approximated using surface temperature *TS* in (K). Calculation of the SWD thus merely requires further conversion of the ZPD, taking into account the incidence angle of acquisitions, in order to obtain the corresponding delay along the LOS. The conversion factor used for the calculation of the theoretical SWD is

$$
\delta^{SW} = \delta^{ZP} \frac{1}{\cos(\theta)}\tag{2}
$$

where *δSW* corresponds to the path delay along the LOS, and *θ* to the incidence angle of SAR observations. The incidence angle of the SAR acquisitions we used in this work (TerraSAR-X descending track 111, see Section 3.3 for details) was inclined 37◦ with respect to zenith path resulting in a conversion factor of 1/*cos*(37) = 1.25.

Finally, the differences between PWV contents estimated for the acquisition dates of the master image and the slave images were calculated (*PWVDIFF* = *PWVmaster* − *PWVslave*), and used to determine the differential SWD (dSWD) that such differential PWV contents theoretically would produce in the gas plumes of each interferogram.

#### *3.5. DInSAR Preparation*

Estimation of the gas plume related radar path delays in interferometric SAR measurements requires several pre-processing steps, involving the removal of large-scale phase contributions related to topography and refractivity changes in the troposphere, the latter of which are commonly referred to as the atmospheric phase screen (APS). First, the raw interferograms were corrected using the baseline information of the interferograms in conjunction with a DEM from the Shuttle Radar Topography Mission (SRTM-1), which provides elevation data at a resolution of one arc second grid spacing (i.e., a relative horizontal precision of ±15 m) and a relative vertical precision of ±6 m. The results are terrain corrected differential interferograms (DInSARs), which however still contain some small-scale topographical errors, due to the coarse resolution of the DEM.

The terrain corrected DInSARs were then unwrapped and further corrected by means of path delay difference maps obtained from numerical weather hindcasts of the WRF, using a horizontal grid spacing of 900 m and a vertical division of 51 eta levels for the simulations. The WRF simulations were initialized with the ERA interim dataset from ECMWF at 06:00 a.m. (UTC) of each SAR acquisition date, and the 900 m fine-resolution domain covering the area of interest further was embedded into two larger coarse-resolution domains using 2700 and 13,500 m grids, respectively, which were used to update the boundary conditions of the inner high resolution domain. This correction was done in order to reduce atmospheric phase contributions, and in particular to avoid altitude correlated distortions known as stratification effect, which is the most prominent atmospheric disturbance in the presence of strong topography [15]. After removal of the linearly stratified atmosphere, the resulting DInSARs, however, still contain phase contributions that were caused by small-scale lateral refractivity variations of a turbulent atmosphere and the volcanic steam plume, which we decomposed as described in the following section.

#### *3.6. DInSAR Decomposition*

To extract the gas plume related phase contribution from our interferograms, we used a wavelet-based DInSAR decomposition (WBDD) technique, which enabled us to integrate our estimates of the PWV contents in the volcanic gas plume into DInSAR analysis and thus allowed for the separate evaluation of gas plume related phase delays and their comparison with respect to all other phase contributions. The WBDD technique [63] was developed to deal with temporally unconnected DInSAR time series subsets and small stack sizes, and does not use the isotropy assumption to mitigate the atmospheric effect. The algorithm is based on Dual-tree complex wavelet transform (*DT*-C*WT*), which is being used in a wide range of applications in the field of signal and image processing [64]. The *DT*-C*WT* calculates the complex transform of a signal using two separate discrete wavelet transform decompositions (tree a and tree b), which may be regarded as equivalent to filtering the input signal with two separate banks of bandpass filters, that separately produce the real and imaginary coefficients of the complex wavelet. The *DT*-C*WT* allows to capture both frequency and location information of a multidimensional signal, i.e., the times and locations at which these frequencies occur.

The WBDD technique requires at least two temporally interconnected interferograms as an input and the workflow is as follows (Figure 3). First, the *DT*-C*WT* representations of all interferograms in the DInSAR time series are derived, i.e., complex wavelet coefficients are determined for each range azimuth position, where the phase of the coefficient determines the exact position of the wavelet in the interferogram and the modulus of the coefficient corresponds to the strength of the wavelet at given position. Second, each of these complex wavelet time series is analyzed using so-called *priors*, respectively *prior time series* containing information on the temporal evolution of any of the processes that possibly contributed to the total phase of the interferograms. That is, we are including prior information on the signal strength variations related to (i) inaccuracies of the DEM, (ii) linear deformation of the ground surface, and (iii) changing refractive properties of the ambient atmosphere and (iv) the volcanic gas plume (details on the determination and compilation of these *priors* are given in the following section). This prior information allows assigning each complex wavelet to their likely causes, depending on their best temporal correlation with the prior knowledge, and therefore enables to decompose the interferograms into signals related to the *priors*. The decomposition thus yields one separate radar propagation path delay estimate for each *prior* that was incorporated in the decomposition analysis of the DInSAR time series.

**Figure 3.** Technical diagram depicting the processing steps of the WBDD algorithm. In a DInSAR time series each DInSAR pixel (respectively range azimuth position) is characterized by a certain temporal evolution of its differential signal strength. Computation of the *DT*-C*WT* representations of each DInSAR enables to capture the temporal evolution of the signal strength of all SAR pixels of the DInSAR time series by means of a small number of complex wavelet coefficients. Similarly each process which causes changes in SAR signal strength can be described by a time series, which reflects the temporal variations of the process (e.g., variations of water vapor contents in the volcanic gas plume, variations in spatial baseline, relative humidity, ground temperature, and pressure). The algorithm uses the time series of these processes as an a priori knowledge of a related possible change in SAR signal strength and assigns the SAR signals to their likely causes by comparison of the *prior time series* with the temporal evolution of SAR signal strength at each range azimuth position, which decomposes the interferograms into different phase screens.

#### *3.7. Priors Used for DInSAR Decomposition*

To model the signal strength variations associated with different sources/processes affecting the interferometric measurements we used the following *prior time series* reflecting the temporal behavior of the corresponding processes. The dSWD time series, which were determined following the descriptions given in Section 3.4, were used as prior information for our algorithm, aiming to determine exclusively those interferometric phase variations, which are characterized by a similar temporal evolution as that of the PWV contents in the gas plume. As a result, we obtained the intersecting set of all gas plume related phase contributions contained in the DInSAR time series. Thus, our gas plume related phase delay estimates exclusively cover the area in which gas plumes of the different interferograms overlap.

The spatial and temporal baseline histories of the interferograms were included as prior information in order to capture small-scale topography related phase contributions, which were caused by errors in the DEM, and respectively by linear ground surface deformation. This works, because the effect of the DEM error typically evolves proportionally to the spatial baseline history of a DInSAR time series [65], and linear ground surface deformation evolves proportionally to the temporal baseline history. The values of the baseline histories were taken as they are and did not require any further computation, since these already represent differential values.

Furthermore, we estimated the refractivity related phase contributions of the ambient atmosphere in order to prevent them from leaking into the gas plume estimate, which involved incorporation of several *priors* in the analysis. In our approach we separately consider non-repeating phase contributions, which are caused by turbulent mixing of the atmosphere, and repeating phase contributions, which are rather related to the vertical stratification of the atmosphere and orographic effects, and thus typically correlate with the underlying topographic relief. *Dirac-function priors* were used to model APS spatial variation at each SAR acquisition, encompassing all non-recurrent phase contributions, which e.g., resulted from the turbulent transport of atmospheric and volcanic water vapor in the lower part of the troposphere, and thus were exclusively present at the times of single SAR observations. For this purpose one *Dirac-function prior* was included for each of the SAR scenes (using positive *Dirac-functions* for master scenes, and negative *Dirac-functions* for slave scenes), providing the benefit that no external information is required for this type of estimate, which in the following will be referred to as so-called *single event APS estimate*. In addition to *single event APS estimation*, we incorporated *prior time series* containing information on the differential temporal variations of relative humidity, surface temperature and surface pressure, which were derived from the detailed WRF simulation (described in Section 3.5) over a single spot on the summit of the volcano. This was done in order to also cover repeatedly occurring small-scale disturbances, which are related to stratification and to recurring formation of localized orographic clouds, and thus had not yet been captured by the *single event APS estimates*. Moreover, this approach decomposes the refractivity related propagation delay into different contributions of the single physical processes, which contribute to the variations of the tropospheric delay (compare with Equation (1) in Smith and Weintraub [66]), and thus enables us to separately examine the small-scale phase delay effects of the different physical processes in the troposphere.

#### **4. Results**

#### *4.1. PWV Contents in the Volcanic Gas Plume of Láscar*

During the days of the Multi-GAS survey in December 2012 and in the investigated period, the activity of Láscar volcano was characterized by quiescent degassing from the active crater (Figure 1). Webcam and visual on-site observations confirmed the frequent wind-blown dispersal of a white gas plume in south-easterly direction over the Puna plateau (Figure 1d). This gas plume therefore very often is also fully captured by the well-located scanning DOAS station. SO2 SCDs measured within the gas plume on the days with TerraSAR-X acquisitions ranged from 2.5 × 1016 to 2.8 × <sup>10</sup><sup>18</sup> molecules·cm−<sup>2</sup> and averaged at 7.2 × 1017 molecules·cm−2. SO2 SCDs exhibited diurnal variations, which were on average slightly enhanced during the early morning and evening hours.

The Multi-GAS measurements conducted during the field survey in December 2012 showed considerable fluctuations of the water contents at nearly constant CO2/S molar ratios, suggesting variable water loss due to condensation of water inside the cooling plume prior to the sampling by the instrument [33]. This in turn resulted in a particularly low average molar H2O/SO2 ratio of 9.37:1, as Multi-GAS instruments exclusively measure gaseous water vapor and do not take into account the LWC in the volcanic cloud. Using the values of Tamburello et al. [33] we thus determined a maximum molar H2O/SO2 ratio of 34:1, by means of which we finally arrived at a water content constituting roughly 90 mol.% of total gas emissions, which is in the range of typical values (85–99 mol.%) for high temperature gases from arc volcanoes elsewhere (e.g., Table 5 in Oppenheimer et al. [67]).

Our PWV estimates were thus obtained by conversion of SO2 SCDs at each scan angle of the plume cross sections using the fixed molar H2O/SO2 ratio of 34:1 (Figure 4a; see Appendix A for details). The amount of water vapor in the centerline of the plume was determined by finding the maximum PWV contents of each DOAS scan, and these were used to constrain the daily average PWV contents for days with TerraSAR-X acquisitions. Estimated daily average PWV contents obtained from the centerline of the plume ranged from 0.007 to 0.015 mm at the times of SAR observations (Figure 4b), translating to SWDs of 0.07 to 0.1 mm (Table 2, columns 3 and 8). Additionally, the average amount of water vapor in each DOAS scan was determined and used to calculate the daily average PWV contents of the bulk plume, i.e., by taking into account the whole plume cross sections at this time. Daily average bulk plume PWV contents ranged from 0.002 to 0.0035 mm (Figure 4b), corresponding to SWDs of 0.02 to 0.03 mm. Differential SWDs from the plume center were used as a *prior time series* (Table 3, column 3) for the estimation of the gas plume related phase contribution in the proximal portion of the gas plume. These values reflect H2O vapor contents at the crater rim, and are thus not representative for PWV contents arriving at the scanning plane of the scanning DOAS (Figure 2a), since they do not take into account downwind evaporation.

Potential evaporation rates were thus determined for the measurement periods of the scanning DOAS (Figure 4c; see Appendix B for details on evaporation calculation), and used to upscale the H2O/SO2 ratio, which was thereby adjusted to increase proportionally to evaporation rates, which in turn strongly depend on wind speeds (Figure 4d) over arid areas. The resulting variable ratio was then used to upscale the SO2 SCDs at each scan angle of the plume cross sections as described above (and in Appendix A), at this time, however, taking into account the effect of downwind evaporation (Figure 4e).

We assume that such adjusted estimates provide a more reasonable approximation of the temporal variations in PWV encountered above various points of the plume cross sections. Estimated daily average PWV contents obtained from the centerline of the plume ranged from 0.8 to 9.6 mm at the times of SAR observations (Figure 4f), translating to SWDs of 6.4 to 77.0 mm (Table 2, columns 4 and 9). Daily average bulk plume PWV contents ranged from 0.2 to 2.5 mm (Figure 4f), which correspond to SWDs of 1.6 to 20 mm (Table 2, columns 5 and 10). The daily average bulk plume PWV contents were used for the determination of the second dSWD *prior time series* (Table 3, column 4), which in turn was used for the estimation of the gas plume related phase contribution in the distal portion of the volcanic gas plume.

**Figure 4.** (**a**) Time series of PWV contents determined for plume cross-sections recorded by scanning DOAS on days with TerraSAR-X acquisitions (18 December 2013 to 16 February 2014). The time scale is not linear. Result without consideration of evaporation; (**b**) Daily average PWV contents in the centerline of the plume (*stippled line*) and the bulk plume (*solid line*), which are representative for the moisture distribution above the crater rim on days with available SAR acquisitions; (**c**) Time series of potential evaporation rates at plume height above the summit of Láscar volcano. Evaporation rates at the times of SAR observations are indicated by *red* (master scenes) and *blue* (slave scenes) *dots*; (**d**) Time series of wind speeds at summit altitude used for calculation of potential evaporation rates; (**e**) Time series of PWV contents in plume cross-sections considering downwind evaporation; (**f**) Daily average PWV content in the centerline of the plume (*stippled line*) and the bulk plume (*solid line*) downwind of the volcano on days with available SAR acquisitions.


*Remote Sens.* **2018** , *10*, 1514

#### *4.2. Gas Plume Related Phase Delays*

PWV contents obtained for the distal portion of the volcanic gas plume (Figure 5a) were used to exemplify how the DInSAR decomposition analysis can be used to map the corresponding refractivity anomalies in the DInSAR time series. To this end, we used the entire DInSAR data set described in Section 3.3, and applied all *priors* listed in Section 3.7. The differential interferograms, which were coarsely corrected according to the steps described in Section 3.5, show a relatively large amount of coherent pixels (Figure 5b). Pronounced phase changes are visible around the summit and on the flanks of the volcano, representing the sum of phase contributions from individual sources including slow and linear ground surface deformation [40], residual topography due to errors in the DEM, and refractivity fluctuations in atmosphere and volcanic gas plume.

The coarsely corrected DInSARs were decomposed (Section 3.6) in order to obtain the intersecting set of repeating gas plume related phase contributions, which are present in all interferograms of the analyzed DInSAR time series. Owing to the different temporal development of PWV contents which were determined for the two different sections of the plume (Figure 4b,f), we accordingly retrieved two contrasting phase delay results for the ascending part of the plume above the crater (Figure 6a) and the downwind portion of the gas plume (Figure 6b). The phase delay map obtained for the *prior time series* of the gas plume related SWDs above the crater rim (Table 3, column 3) displays a weak altitude correlated shortening of the delay over an area which is largely confined to the summit region of the volcano (Figure 6a). The phase delay map obtained for the downwind portion of the gas plume (Figure 6b) in contrast shows a lenticular pattern that indicates lengthening of the radar propagation path over an area that agrees well with the most common eastward plume transport directions observed during spring and late summer of that period (see Appendix A, Figure A1c,d).

DInSAR decomposition further yielded phase delay maps for each of the other *priors* that were included in the analysis. A detailed description of these maps (which are displayed in Figures 6 and A4) would however go beyond the scope of this paper and thus can be found in the Appendixs E and F.

**Figure 5.** (**a**) PWV contents in plume cross-sections obtained for the downwind portion of the volcanic gas plume on acquisition dates of SAR master and slave scenes, and (**b**–**e**) corresponding DInSAR maps (10 × 10 km in azimuth and range direction). Scales of the DInSAR maps indicate range change in units of mm, and are unique to each image using the same scale bounds for the sake of comparability with other phase contributions. Incoherent areas are masked out; (**b**) Coarsely corrected DInSAR maps used for decomposition analysis and retrieval of the gas plume estimate. DEM and major APS contributions derived from WRF were removed; (**c**) Phase difference maps of the gas plume estimates obtained for the downwind portion of the volcanic gas plume. Corresponding theoretical dSWDs are indicated in the upper right corner of each map; (**d**) DInSAR maps from 5b, where the phase contributions of the gas plume were removed; (**e**) Linear deformation estimates obtained from the refined correction using the delay estimates obtained from WBDD analysis.

**Figure 6.** (**a**–**h**) Delay correlation maps depicting the estimated repeating patterns of SAR delays obtained from *priors* included in the WBDD analysis. Location of active crater is indicated by the *red circle*; (**a**) Delay estimate for the gas plume aloft the crater (**b**–**d**) Delay estimates for the downwind portion of the gas plume (**b**) with and (**c**) without additional removal of pressure and temperature dependent phase screens by including/excluding respective *priors* in the WBDD analysis; (**d**) Gas plume estimate as in c), however omitting the SAR observation of 3 January 2014. Scales correspond to estimated delay (mm) per theoretical delay (mm). Scale bar limits of the gas plume estimate obtained for the proximal part of the plume above the crater indicate that theoretical delays are by a factor 100 smaller than the delay estimate. Upper and lower limits of the scale bars of the estimates obtained for the downwind portion of the gas plume are equal to unity, indicating concurrence of theoretical and estimated delay. The red to yellow colored signature in the lower right corner of the image indicates a lengthening of the delay, where the effect of H2O emissions is large. The affected area corresponds to the most common plume transport directions (see Appendix D, Figure A3a,b) and locations where the volcanic plume regularly touches the ground (fumigation); (**e**) Estimates of surface temperature and (**f**) surface pressure related delays. Scales indicate mm estimated delay per Kelvin, respectively hPa of the input *prior*; (**g**) Delay estimate obtained for the relative humidity *prior* (**h**) Delay estimate for the spatial baseline *prior*; (**i**) Scatter plot of theoretical delay (mm) versus determined signal strength of the gas plume contribution in the DInSARs used for the gas plume estimate in (**b**). Each asterisk symbol represents one of the DInSARs and corresponding slave dates (mm-dd) are indicated.

#### *4.3. Isolation of the Gas Plume Related Signal in DInSARs*

The phase delay map obtained for the downwind portion of the gas plume (Figure 6b) was used to determine the gas plume related phase contributions in each DInSAR image of the time series. For this purpose the theoretical dSWDs, which were used as a *prior* for this estimate (Table 3, column 4), were multiplied by the respective phase delay estimate, which resulted in simple phase difference maps depicting the repeating gas plume related phase contribution present in each of the DInSAR maps (Figure 5c). The phase difference maps use the same scale bounds as those used in Figure 5b in order to display the relative effect of the gas plume related phase delay on the interferometric measurements with respect to other phase contributions. Each of these maps reveals a lenticular phase difference pattern, which equals the pattern that was already observed in the corresponding phase delay map (Figure 6b), and which persists throughout the entire DInSAR time series. As a consequence of the mathematical operation used to retrieve these maps, the amplitude of the pattern varies in proportion to the a-priori constraints on the differences in gas plume PWV content and is most pronounced in those DInSARs, in which the differential PWV contents determined for the epochs of master and slave SAR acquisitions were the highest (compare Figure 5a,c). Moreover, it is interesting to note that the contrast between plume related and atmospheric signals was more pronounced in the DInSAR images of *subset 01*.

As a next step, the gas plume related phase difference maps depicted in Figure 5c were subtracted from the coarsely corrected DInSAR maps depicted in Figure 5b in order to mitigate the gas plume related phase contribution. The resulting DInSAR maps (Figure 5d) differ clearly from the coarsely corrected DInSAR maps, which still contain the gas plume related signal. As expected, the removal of the gas plume related signal is most obvious in the DInSARs of *subset 01*, where the mitigation in two cases (DInSARs 1 and 3) even resulted in a change to the opposite delay direction over the area which was affected by the gas plume (Figure 5d, rows 1 and 3).

#### *4.4. Residuals of the Refractivity Related Phase Delay and Ground Deformation*

To verify the efficiency of our decomposition analysis, we now examine the residual phase, which persists when all estimated refractivity and DEM error related phase contributions have been removed from the DInSAR maps. To this end, the complete set of refractivity and DEM error related phase delay estimates obtained from WBDD analysis as shown in Figures 6 and A4 (and further detailed in Sections 4.2 and 5.3 and the Appendix) were multiplied with the values of their respective input *prior time series* to obtain phase difference maps corresponding to each of the DInSARs. Similarly the *single event APS estimates* of individual SAR acquisitions (Figure A4) were combined to form phase difference maps for each interferogram by subtraction of the slave APS from the master APS. The resulting phase difference maps were then subtracted from the coarsely corrected DInSARs (Figure 5b) obtaining interferograms, which mainly contain the residual phase and topographic contributions attributed to a deforming ground surface (Figure 5e). This refined correction more or less efficiently removed small-scale phase delay patterns, which were caused by the volcanic gas plume and a turbulent atmosphere, as well as small-scale DEM errors related to acquisition geometry, respectively spatial baseline, leaving only some faint and blurry signatures of low amplitude.

In between the blurry residuals of non-deformation related phase contributions a putative small-scale ground deformation signal can be observed in the central of the three north-eastern summit craters of Láscar. The signature occurs in each of the DInSARs, and indicates a progressively subsiding crater floor confined to an arcuate fault, where Pavez et al. [39] detected co-eruptive subsidence of the crater floor during the 1995 eruption, and where Richter et al. [40] identified linear subsidence during the period 2012–2016.

#### *4.5. Validation of the Gas Plume Signal Estimate*

The phase delay estimates which were obtained by decomposition analysis can be validated by means of relating the position and strength of the estimated DInSAR signal to the input prior information. We thus compared the signal strength of the lens-shaped gas plume pattern in each of the DInSAR maps (Figure 5b) to the strength of the theoretical dSWDs (Table 3, column 4), which served as a *prior* for the corresponding gas plume estimate (Figure 6b). For this purpose the signal strength *S* of the gas plume related phase contribution was determined for each interferogram using the amplitude of the gas plume pattern in the phase delay estimate obtained for the downwind portion of the plume (orange color in Figure 6b) by means of following Equation (3).

$$S = \operatorname{sgn}(\overrightarrow{X} \cdot \overrightarrow{Y}) \sqrt{|\overrightarrow{X} \cdot \overrightarrow{Y}|} \tag{3}$$

here *sgn* denotes the sign function and <sup>→</sup> *X*· → *<sup>Y</sup>* <sup>=</sup> <sup>∑</sup>*<sup>i</sup> XiYi* is the dot product of vectors <sup>→</sup> *<sup>X</sup>* and <sup>→</sup> *<sup>Y</sup>*, where <sup>→</sup> *X* corresponds to the amplitude of the gas plume estimate obtained from WBDD and <sup>→</sup> *Y* is the associated measured signal in each of the DInSARs, which were used in the analysis. The comparison revealed a moderate positive linear relationship between signal strength *S* and theoretical strength (Table 3, column 4), which is fairly well represented by the regression line depicted in Figure 6i, since 61% (R2 = 0.61) of the total variation in measured signal strength can be explained by the linear relationship between theoretical strength and signal strength. Thus, in principle, this regression line can be used to determine the gas plume related phase delays even of those DInSARs for which no corresponding ground-based gas emission measurements are available.

#### **5. Discussion**
