*5.1. Why the Distinction between Periods of Volcano Deformation and Enhanced Degassing is Important*

This work constitutes the first serious attempt to quantify gas plume related phase delays in differential interferometric radar measurements of a persistently degassing and deforming volcano. We propose a method that enables us to attain gas plume detection by means of satellite-based radar interferometry and demonstrate that gas plume related phase delays in DInSARs can be isolated and mapped, if a-priori constraints on the gas plume PWV contents obtained from independent gas measurements are included in the analysis of DInSAR data.

The distinction between periods of local volcano deformation and enhanced degassing can be of vital importance both for early warning purposes and for a better understanding of volcano-tectonic processes. Degassing and deforming volcanoes often display a variety of deformation sources, which reflect the complex interplay between (i) magma emplacement/retreat, (ii) accumulation/discharge of volcanic gases and (iii) changes in hydrothermal activity, as e.g., was observed at Campi Flegrei in Italy [68,69], or at Lastarria in Northern Chile [70]. Retention and loss of magmatic and hydrothermal volatiles play a central role in virtually any of the processes that lead to volcano deformation. Gas emission measurements on deforming volcanoes thus provide essential auxiliary information on the ultimate cause of individual deformation observations [71,72], and deformation measurements likewise may deliver explanations for observed changes in the amount and composition of gas emissions [73].

#### *5.2. How Methodological Limitations Can be Turned into Benefit*

Area-wide precise ground deformation measurements by means of InSAR require adequate mitigation of atmospheric phase contributions, which may be challenging in small SAR datasets. In particular, repeatedly occurring spatially correlated errors, which e.g., may be caused by orographic clouds and volcanic gas plumes, are difficult to remove [74,75]. Insufficient resolution of meteorological and DEMs used for correction of DInSAR data adds further uncertainties to retrievals of ground

deformation. Moreover, most modern SAR satellites revolve Earth on a sun-synchronous dusk-dawn orbit, and record their data along the day-night boundary, i.e., at a time, when emission rates of continuously degassing volcanoes typically are the highest [76–78], which in turn increases the risk of having disturbances due to volcanic water vapor emissions in InSAR measurements. Such noise may, however, be transformed into data, as soon as phase contributions of different sources can be distinguished, which can be accomplished by integration of a-priori constraints.

In our approach we derive prior constraints of volcanic H2O emissions using measurement techniques, which are commonly used in volcano monitoring, hence providing the benefit of data availability for a comparatively large number of volcanoes. Volcano monitoring generally demands rapid response, in order to assess the volcanic hazard in a timely manner, and the specific information which would be relevant to identify more rapid ground deformations often may be confined to a limited number of available SAR observations. Therefore it seems highly desirable to develop a method that allows removing disturbances of any kind even from small sets of interferograms, where conventional methods as persistent scatterer interferometry and the small baseline subset algorithm would fail. The presented method addresses the problems of the different measurement methods mentioned above, and exploits the fact that DInSAR analysis allows for mapping of water vapor distribution.

In the following subsections, we will discuss the phase delay effects that can be observed in the phase delay maps which were obtained for the proximal and distal sections of the plume, and relate our observations to variations in gas emission rates. Furthermore, we will examine the relative and absolute effects of these phase delays with respect to (seasonal) variations in water vapor contents of the ambient atmosphere. Finally, we will assess the benefits and drawbacks associated with our approach.

#### *5.3. Phase Delay Effects above the Crater and Downwind of Láscar*

We now elucidate the contrasting phase delay effects that are observed in the phase delay estimates obtained for the proximal and distal portions of the gas plume (Figure 6a,b).

The phase delay map determined for the downwind portion of the gas plume (Figure 6b) shows a lenticular pattern that indicates lengthening of the radar propagation path due to enhanced moisture 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 D, Figure A1c,d). Westward drifting plumes, which generally may have occurred during early summer (see Appendix D, Figure A1c,d and Figure A2), did obviously not occur at times of the SAR acquisitions we used to obtain our gas plume related phase delay estimates, as these would have prohibited detection of a repeatedly occurring phase delay pattern. Moreover, such westward drifting plumes would less likely leave any repeating traces in the SAR signal, due to a less pronounced humidity/refractivity contrast between plume and surrounding atmosphere in early summer. This is further supported by the fact that easterly winds over the Atacama region typically have a lower velocity (see Appendix D, Figure A2), and tend to be less stable with respect to their direction (see Appendix D, Figure A1c,d), resulting in unsteady plume transport directions, which additionally impedes formation of pronounced repeating patterns. Wind speeds were generally larger during periods with eastward plume transport (spring and late summer; see Appendix D, Figures A1a and A2), causing plume transport to be more turbulent, and thus the volcanic plume was very likely less condensed, and instead contained more vapor during that period.

Additional phase delay estimates of the distal portion of the gas plume were processed without (Figure 6c,d) additionally including prior information on the temporal behavior of surface temperature and pressure in order to show the effect of APS mitigation by the corresponding temperature and pressure phase screens (Figure 6e,f). For the computation of the estimate displayed in Figure 6d we furthermore reduced the number of input DInSARs omitting the interferogram that was formed from SAR acquisitions of 18 October 2013 and 03 January 2014 (Figure 5, row 5). Prior information on relative humidity variations were included in all of the mentioned cases, yielding a relative humidity related phase delay estimate (Figure 6g), which was mitigated from all other phase delay estimates. Comparison of the resulting gas plume phase delay estimates clearly shows that the repeating patterns of the gas plume signal (Figure 6c,d) can readily be reduced (Figure 6b), when a surface temperature and air pressure correction (Figure 6e,f) is included in the WBDD analysis. The gas plume related phase delay estimate depicted in Figure 6b thus clearly represents the result in which the separation of gas plume related phase delays from other phase contributions was most successful.

The phase delay map obtained for the prior time series of the gas plume related SWDs above the crater rim (Table 3, column 3) shows 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 location of this phase delay signature is thus consistent with the assumption that the dSWDs, which were used as a *prior* for this estimate (Table 3, column 3), are representative for the temporal evolution of SAR signal strength variations measured above the crater. The observed shortening of the delay, however, contradicts the hypothesis that this signature can be ascribed to enhanced moisture over the respective area, since a water vapor field would produce a lengthening of the delay (compare Equation (1) in Smith and Weintraub [66]). Furthermore, the phase delay patterns in this estimate strongly resemble the patterns observed in the DEM error related phase delay estimate (Figure 6h), which was obtained using the temporal history of the spatial baseline (Table 3, column 8), suggesting that the phase delays, that the algorithm attributed to the gas plume above the crater, are in fact rather related to errors in the DEM (detailed in Appendix F). This suspicion is further substantiated, if one compares the shape of the two corresponding *priors*, which strongly resemble each other in appearance. As a consequence of this similarity, the signal obviously was assigned to both *priors* by the WBDD algorithm, since this signal does not perfectly match to any of the *priors*. This phase delay estimate thus is a good example for a case in which the algorithm failed to isolate possible gas plume related phase delays from other phase contributions. Thus, care must be taken to avoid misinterpretation of model results.

#### *5.4. Gas Emission Rates from Láscar Volcano*

Gas emission rates from Láscar volcano were determined and compared to our PWV estimates (Figure 4b,f), in order to illustrate how much water vapor has been emitted from the volcano to produce the corresponding PWV contents in the gas plume, and to relate our observations to variations in volcanic degassing activity. SO2 emission rates were calculated from scanning DOAS measurements using wind field information obtained from GDAS1 soundings provided in the web-based archive of the National Oceanic and Atmospheric Administration's (NOAA) Air Resources Laboratory [79]. These were then upscaled by the H2O/SO2 molar ratio of 34:1 obtained from Multi-GAS, in order to retrieve the corresponding H2O-fluxes (Figure 7a).

Assuming a constant H2O/SO2 ratio and a constantly strong condensation of the emitted water vapor, the H2O (SO2) emission rates ranged from 5 to 100 kt·day−<sup>1</sup> (150 to 2900 t·day<sup>−</sup>1) and averaged at 12 kt·day−<sup>1</sup> (350 t·day−1) during the considered period. Based on these assumptions we estimate that the total errors of our upscaled H2O emission rates likely were around −15 to +45% and thus even more skewed towards underestimation than the error envelopes of corresponding SO2-fluxes, which were constrained to about −10 to +30% (see e.g., [80] for error determination). Gas emission rates of October and November 2013 showed more pronounced variations and were on average slightly enhanced (about 19 kilotons of H2O and respectively 550 tons of SO2 per day) with respect to the rest of the period (10 kilotons of H2O and respectively 300 tons of SO2 per day). The enhanced degassing activity of October and November further was accompanied by elevated heat dissipation from the main active crater as was suggested by intermittent night-time webcam observations of illuminated gas plumes (obtained on several days during the first half of October and lastly on 20 November 2013; see Figure 7a), which indicated the presence of incandescent material inside the active crater [81]. Such observations support the repeated influx of fresh magma from depth feeding the high-temperature gas exhalations at the surface of the crater bottom, thus suggesting a concurrent increase of magmatic gas emissions and progressive drying of the ambient hydrothermal system (e.g., [41]), which may have altered gas ratios temporarily. Contemporaneous measurements of SO2

emission rates confirmed this suspicion and appropriately showed sharp increases over the course of each of these days (Figure 7a).

**Figure 7.** (**a**) H2O emission rates from Láscar volcano along with a sequence of Aster, Landsat-8 and EO-1 ALI scenes depicting snow and cloud coverage (upper row: false-color composites of SWIR, NIR and visible red spectral bands; lower row: true-color images using a combination of visible red, green and blue spectral bands). Days with incandescence observations are indicated by *stippled grey vertical lines*. Precipitation events that were recorded in Toconao (on 10 November 2013, 17 January 2014, and 26 January 2014) are indicated by *vertical blue lines*, including information on amount and duration of the events; (**b**) Background atmospheric PWV contents estimated from GDAS1 soundings above Láscar volcano (*red curve*), compared to radiometer measurements conducted at APEX (*blue curve*). Good agreement between the two curves reflects similar weather conditions caused by a similar morphologic exposition of both sites, and may additionally be attributed to the coarse spatial resolution of the GDAS1 data. PWV estimates for the times of SAR observations are indicated by *red* (master scene) and *blue* (slave scene) *dots*, respectively.

Bearing these observations in mind, we compared the gas emission rates with the estimated PWV contents of the gas plume (Figure 4b,f) revealing that gas emission rates evolved differently than the PWV contents obtained for the gas plume above the crater (Figure 4b), but displayed a very similar temporal evolution as that of the PWV contents, which were obtained for the downwind portion of the gas plume (Figure 4f). The latter thus likewise were enhanced during the period of increased volcanic activity.

#### *5.5. Estimation of Background Atmospheric PWV Contents*

Atmospheric disturbances due to volcanic gas plumes are intuitively expected to be more pronounced during periods which are characterized by strong variations of gas emission rates and low ambient atmospheric humidity. Determination of PWV contents in the local atmosphere was thus done to be able to approach the relative and absolute effects of the plume on the specific SAR acquisitions. In the Atacama region the driest conditions are typically found during the period late March to early December, and between local midnight and 11:00 a.m. [82]. Background PWV contents in the atmosphere above the volcano were retrieved from vertical atmospheric GDAS1 profiles provided by NOAA (see Appendix C for equations). These estimates were also cross-checked with PWV estimates obtained from radiometer measurements (Figure 7b) provided by the Atacama Pathfinder EXperiment (APEX), which is located about 30 km north of Láscar at a very similar altitude of 5100 m above sea level (m.a.s.l.), thus reflecting similar weather conditions.

Because PWV is a measure of how much water is available for potential rainfall, it dictates precipitation intensity, which thus reversely can be used to infer accompanying PWV. To give an impression of when to expect enhanced PWV contents and to further complete the picture of the hydrological cycle of the study area, we will first briefly describe the characteristic seasonal patterns according to which precipitation typically occurs at Láscar, and narrow down the periods during which precipitation actually was encountered in the course of the observation interval of our case-study.

The considered period covers parts of the "dry season", which is characterized by dry and cold prevailing westerly synoptic winds, and is spanning most of the year (here October, November and second half of February are included), as well as the "wet season" ranging from December to February, which is characterized by less stable atmospheric conditions, resulting in more humidity being transported to the Andes by frequently recurrent warm easterly winds [83]. As a consequence of this easterly moisture source, precipitation on the leeward slopes of the Puna plateau exhibits a marked seasonality in particular at elevations below 3500 m.a.s.l., where about 90% of the annual precipitation (<20 mm·year−1) falls during the (austral) summer months, that is from December to March [84]. The associated rainfall events typically are very short, but can become quite intense. During the considered period, a weather station (operated by the Instituto de Investigaciónes Agropecuarias) located at 2500 m.a.s.l. in Toconao, a small village roughly 30 km NW of Láscar, recorded 4 events of precipitation which amounted to a total of 5.7 mm (equivalent to liters per square meter) within only 3.6 h (vertical blue lines in Figure 7a). At elevations above 3500 m.a.s.l. the situation is slightly different and precipitation (<200 mm·year−1) occurs even during the dry winter period, mostly in the form of snow with the greatest snowfall frequency happening to occur near the Tropic of Capricorn [85], where also Láscar coincidentally is located. Láscar volcano thus commonly is almost completely (except for the hot bottom of the currently active crater) covered with several meters of snow during the "dry season", which however melt and sublimate due to the warm winds of the "wet season". Such snowmelt events hence constitute an important source for groundwater recharge, which potentially affect the hydrothermal activity of Láscar. The temporal development of the snow and cloud cover encountered at Láscar during the period of our case-study is illustrated by a sequence of satellite images (Aster, Landsat-8, and EO-1 ALI) in the top panel of Figure 7a in order to give a visual impression of the situation associated with different epochs of the temporal development of atmospheric moisture.

Measured and estimated atmospheric PWV contents were at average below 1 mm during October and November 2013, sharply increased to 5 mm in December 2013, stayed enhanced through the first half of February 2014, and finally decreased again to 1 mm in the second half of February 2014 (see Figure 7b and also Table 2, column 6). The comparatively low PWV values determined for the times of SAR observations used in this work reflect the pronounced atmospheric transparency around sunrise, and the resulting SWDs (Table 2, column 11) are in good agreement (R2 = 0.96; see Appendix G for detailed statistics) with the amplitudes of associated phase delays observed in the corresponding DInSARs (Figures 5b and A4).

It is interesting to note that the apparent SO2 (H2O) emission rates generally were particularly low whenever the ambient atmospheric water vapor levels were elevated, which is clearly visible especially during the "wet season" (compare Figure 7a,b). This observation suggests the presence of low hanging clouds beneath the volcanic gas plume, which may have caused attenuation of the absorption signal due to contributions of scattered light that did not pass through the plume [86]. Moreover, a fraction of the degassing SO2 may have been scrubbed by dissolution into the aqueous phase e.g., during possible events of fresh water influx into the hydrothermal system [87], though this would most likely only have a marginal effect on the emissions of the hot degassing dome of Láscar, which is less probably being deeply infiltrated by permeating waters that interact with ascending magmatic fluids, than the peripheral low-temperature vent sites. Another mechanism, which may be considered to explain the observed diminution of SO2-fluxes is the dissolution of emitted SO2 in the water droplets of low hanging clouds [88], which, however, usually is not particularly effective over a short distance, and the associated scavenging rate of SO2 emissions is generally estimated being on the order of merely a few percent of total emitted SO2 per hour. This in turn leads to relatively long atmospheric lifetimes of SO2 emissions, which were shown to range between several hours to days according to weather conditions [89].

Estimated PWV contents inside the volcanic gas plume (Figure 4b,f) and modeled atmospheric PWV contents outside the plume (Figure 7b) were compared in order to assess the humidity/refractivity contrast between ambient atmosphere and volcanic plume. The comparison revealed that PWV contents at the downwind plume centerline (0.8 to 9.6 mm) were generally larger than the PWV contents of the ambient air (<1 mm) during periods with dry atmospheric conditions, i.e., most of the year. During early summer the situation was different, and the downwind plume PWV content (<1 mm) was less contrasting with respect to the PWV content of the atmospheric background (2–5 mm). PWV contents of the downwind bulk volcanic plume (0.2 to 2.5 mm) were accordingly much less contrasting most of the year when compared to background atmosphere, and were slightly higher than ambient atmospheric humidity only during the first three SAR observations of the time series (on 18 October 2013, 20 November 2013, and 01 December 2013). Future studies exploiting a larger InSAR dataset of Láscar might therefore well capture these seasonalities in gas plume and atmospheric PWV contents, the latter of which were already identified at Láscar [39] but also elsewhere, e.g., at Campi Flegrei and Vesuvius [90]. PWV contents determined for the proximal part of the plume above the crater rim were in general negligible with respect to PWV contents of the ambient atmosphere (0.007 to 0.015 mm on the centerline of the plume, and 0.002 to 0.0035 mm averaged over the bulk plume).

The corresponding SWDs obtained for the proximal part of the plume above the crater rim ranged from 0.07 to 0.1 mm in the centerline of the plume (Table 2, column 8) and from 0.02 to 0.03 mm in the bulk plume and thus were negligible compared to the InSAR accuracy. The SWDs in the more distal part of the plume downwind of the crater, however, were on the order of background atmospheric SWDs (5.3 to 41.5 mm; compare with Table 2, column 11), and ranged from 1.6 to 20 mm in the bulk volcanic plume and from 6.4 to 77.0 mm along the plume centerline (Table 2, columns 9 and 10), which is more than large enough to produce a detectable effect in interferometric measurements, as can be seen in the corresponding phase delay maps (Figure 5c).

Our results thus collectively suggest that refractivity changes in volcanic gas plumes may, particularly under dry climatic conditions, have a significant effect on DInSAR measurements, hence

requiring extreme caution, when alleged deformation signals are detected above or downwind of visibly degassing volcanoes, especially if these are located at high elevations.

#### *5.6. Limitations of the Proposed Method*

Our approach was hampered by several practical limitations, which were mainly related to spatial data coverage and temporal coincidence of the different measurement techniques that we combined.

Simultaneous measurements of Multi-GAS and DOAS instruments were unfortunately not available for the period considered for DInSAR analysis, and thus our PWV estimates had to rely on the gas concentration measurements obtained during a short Multi-GAS survey conducted about one year prior to the period considered here. This is of particular relevance, since Láscar exhibited elevated SO2 emission rates and incandescence during the first two months of this period, suggesting a slightly enhanced magmatic activity, which likely was accompanied by temporarily altered gas ratios. At the same time, the snow cover on the summit of the volcano was observed to be reducing strongly (satellite images in Figure 7), which may be an indication for concurrent fresh water influx into the peripheral hydrothermal system. Thus, there is evidence for both, elevated emissions of magmatic SO2 and increased hydrothermal H2O contributions of meteoric origin. The assumption that molar H2O/SO2 ratios stayed more or less constant throughout the entire period of our case study therefore may not be perfectly appropriate. Availability of data from a stationary Multi-GAS instrument certainly would have helped to better characterize the temporal variations of molar H2O/SO2 ratios, and thus to better quantify the water vapor contents at the instance of gas sampling by the Multi-GAS.

Moreover, volcanic gas plumes typically contain variable amounts of condensed liquid water in addition to gaseous water, which are not taken into account by the measurements of Multi-GAS instruments, and thus unequivocally lead to an underestimation of the water vapor contents, which are expected to increase in the downwind portion of the volcanic gas plume due to evaporation. Such an underestimation could be avoided, if water vapor contents were mapped in situ across the whole extent of the gas plume, e.g., by means of airborne gas composition measurements, though this approach would not be very practicable but rather labor and cost-intensive in the view of satellite revisit times of 11 days.

Additionally, the DOAS measurements are generally restricted to daylight conditions and thus deviate from SAR acquisition times by up to 2 h (Figure 2b). Due to this deviation in acquisition times and because the proposed DInSAR decomposition analysis requires only one representative PWV value for each SAR image (see Sections 3.6 and 3.7), we thus used daily average gas plume PWV contents, in order to determine the SWD that such PWV contents theoretically would produce in the LOS of a SAR image. In the view of strong diurnal degassing variations, which are characterized by pronounced gas emission maxima around sunrise and sunset [77,78], such time averaged values however inevitably lead to systematic underestimation of the water content in the volcanic plume at the times of the SAR observations.

Furthermore, our analyses have shown that different locations within the gas plume may exhibit a different temporal evolution of the PWV contents, which may appropriately be captured by inclusion of corresponding *priors*. However, we point out that modeling the evolution of water vapor contents over the entire extent of the gas plume even by means of two *prior time series* still is a simplification. Prudent sampling of the gas emission data in time and space hence constitutes an important prerequisite in order to obtain a-priori information, which yields feasible results from DInSAR decomposition analysis.

This becomes particularly clear, if one notes that the position of the plume center varied during the measurement period, because the gas plume of Láscar somewhat meanders within the scan range of the scanning DOAS. Transport of the plume changed between north-easterly and south-easterly directions approximately every 10 min, which roughly corresponds to the duration of 1 complete DOAS scan. This resulted in different plume positions in consecutive scans, which is why the sequences of consecutive scans in Figure 4a,e have a corrugated appearance. We thus assume that the resulting

daily average plume center statistic may not be representative for the effect that the bulk plume has on SAR measurements further downwind of the volcano, since the gas distribution along plume cross sections is highly irregular when the plume meanders. For this reason, we chose to use the PWV contents of the bulk gas plume instead of those from the plume center in order to estimate related phase delay in the downwind portion of the gas plume.

Further limitations of our approach lie within the applicability of the WBDD technique. The method makes extensive use of ancillary data, which have to be made available to obtain reasonable results. To successfully isolate the plume-induced signal from total phase, it is necessary to additionally include the *priors* of all potentially interfering signals in the analysis, i.e., the temporal and spatial baseline histories, which were used as *priors* for the estimation of phase delays related to DEM error and linear deformation, as well as several *priors* defining the atmospheric phase delay of each acquisition. We consider that, even when all required *priors* are being used, the estimated gas plume related signal still may additionally contain phase information related to any other interferometric signal that coincidently correlates with the gas plume related signal (as was shown in Section 5.3 for the phase delay estimates depicted in Figure 6a,h). In the presence of a non-linearly deforming ground surface it may additionally be necessary to include prior information on the temporal behavior of the associated deformation signal, which e.g., could be derived from deformation measurements of independent ground-based methods, such as GPS, or tiltmeters. Last, but not least, the DInSAR decomposition analysis exclusively yields reasonable results for the gas plume estimate when there is significant overlap between gas plumes of different SAR acquisitions, since "solitary" plumes would be captured by the *single event APS estimates* causing a gap in the temporal evolution of the repeating pattern. To capture the whole gas plume at the times of SAR acquisitions instead of only the intersecting area of multiple different gas plumes definitely would have been an asset.

#### *5.7. Advantages of the Proposed Method*

The WBDD technique provides a framework for data fusion, enabling to integrate a variety of external information, such as weather data, results from gas emission monitoring or even independent deformation measurements into the analysis of DInSAR time series. These time series may be very short and it is shown that the combination of only three SAR acquisitions already enables us to generate high-quality DInSARs including estimates of their corresponding phase screens (as demonstrated in Appendix H). The algorithm does not require any information on the location of the plume, in order to find the related "disturbance" in the interferometric signal. This is of particular advantage to the application described here, since we lack precise information on the heading of the plume at the times of SAR observations, apart from a rough approximation that can be derived from webcam images and wind directions offered by the weather models. Also, no information about the interferometric signal strength, rather merely the shape of the prior time series is used, because prior information about the relations between signal strength and amplitudes of the prior time series is allotted by chance, which is why the prior information can be scaled arbitrarily without influencing the location of the associated filtered phase delay patterns. Correspondingly, the radar path delay estimates can be validated by the position and strength of the estimated signal in relation to the corresponding prior information (Section 4.5). The gas plume is transported by advection, therefore has a direction, and is correspondingly anisotropic. For this reason, the WBDD technique arguably is the best choice for the detection of gas plumes in interferometric measurements. Additionally, the ability to use temporally unconnected time series subsets enables to use only those DInSARs which have a very small spatial baseline. This helps avoid unwrapping errors due to erroneous DEMs.

#### **6. Summary and Conclusions**

Here, we presented a novel time-space-based filtering method that allows for the isolation and mitigation of virtually any conceivable non-intrinsic error source that affects interferometric ground deformation measurements, by means of including prior constraints on their temporal behavior in the DInSAR analysis. Including *prior* information on gas plume related refractivity variations furthermore enabled us to identify and map spatially correlated repeating interferometric phase delay patterns caused by refractivity changes within volcanic gas plumes that were drifted in a similar direction during several subsequent SAR acquisitions. Apart from potentially being a new tool to detect volcanic gas plumes in DInSAR data, the WBDD algorithm thus provides the possibility to mitigate the plume-induced phase delay, where deformation measurements are the main purpose of monitoring.

Similarly to the small baseline subset algorithm our technique operates on a pixel-by-pixel basis on areas that show sufficient coherence throughout the entire DInSAR time series. The algorithm requires only a relatively short time series comprising at least 2 temporally interconnected DInSARs, respectively 3 observations in order to capture the plume related phase delay, and when more than 2 interferograms are available, these do not necessarily have to be temporally connected. Moreover, the method allows for iterative determination of PWV contents in the volcanic gas plume by matching the estimated phase delays with theoretical phase delays derived from ground-based gas emission measurements, which resulted in reasonably realistic values for the gas plume of Láscar volcano.

To our knowledge this is the first time that phase delay effects in volcanic water vapor emissions were quantitatively investigated by means of radar interferometry. Examination of PWV contents and associated phase delay effects in the volcanic gas plume of Láscar yielded daily average bulk plume PWV contents of 0.2 to 2.5 mm water column, which would generate plume wide excess path delays in the range of 1.6 to 20 mm. The corresponding H2O emission rates ranged from 5 to 100 kt·day−<sup>1</sup> and remarkably displayed a similar temporal behavior as the PWV contents that were obtained for the downwind portion of the gas plume. Our observations consider that the locations of the interferometric patterns in the phase delay estimates are consistent with determined locations of the volcanic plume, and variations in phase delay amplitudes are in good agreement with variations in the strength of the emission source and humidity of the ambient atmosphere. Therefore there is ample evidence, that the DInSAR phase variations detected by the WBDD algorithm are indeed related to refractivity variations in the volcanic plume. Additional independent data on plume location and chemistry at the time of SAR acquisitions would help to further corroborate our results as a proof of concept.

We have illustrated that integrating a-priori information from gas emission measurements and meteorological data into the analysis of DInSAR time series provides previously unknown possibilities to efficiently decompose the radar signal into different phase contributions. The implications of this new technique are wide, ranging from improved deformation measurements to simultaneous quantification of volcanic outgassing activity, DEM error estimation, and meteorological applications, such as mapping and quantification of turbulently transported atmospheric PWV in each SAR acquisition. Our results furthermore encourage the use of data from ground-based gas emission monitoring networks, such as NOVAC for calibration and validation of satellite-based measurements on active volcanoes.

**Author Contributions:** Conceptualization, S.B. and F.-G.U.; Methodology, S.B. and F.-G.U.; Software, F.-G.U.; Validation, F.-G.U. and S.B.; Formal Analysis, F.-G.U. and S.B.; Investigation, S.B. and F.-G.U.; Data Curation, S.B. and F.-G.U.; Writing-Original Draft Preparation, S.B., F.-G.U., T.H.H. and T.R.W.; Writing-Review & Editing, S.B., T.H.H., T.R.W., and F.-G.U.; Visualization, S.B. and F.-G.U.; Supervision, T.R.W. and T.H.H.; Project Administration, T.R.W. and T.H.H.; Funding Acquisition, T.R.W. and T.H.H.

**Funding:** This project was funded by the Helmholtz Association through the "Remote Sensing and Earth System Dynamics Alliance", which was realized in the framework of the Initiative and Networking Fund, and by the GEOMAR Helmholtz-Centre for Ocean Research Kiel (project-id: HA-310/IV010). This is a contribution to VOLCAPSE, a research project funded by the European Research Council under the European Union's H2020 Programme/ERC consolidator grant no. ERC-CoG 646858.

**Acknowledgments:** The authors would like to thank Claudia Bucarey from Observatorio Volcanológico De los Andes del Sur (OVDAS) for installation and maintenance of the permanent scanning DOAS instrument at Láscar volcano, providing an invaluable gas emission data set. TerraSAR-X High Resolution Spotlight scenes were provided by DLR and were tasked and acquired by the proposal WA1642. Further acknowledged is the Atacama Pathfinder EXperiment for providing archived weather and radiometer data, and the Instituto de Investigaciónes Agropecuarias for providing the precipitation data of their weather station in Toconao on the server of the Dirección Meteorológica de Chile (Meteochile). The Aster, Landsat-8 and EO-1 ALI images used in this study are made available by the U.S. Geological Survey through the USGS EarthExplorer platform (https://earthexplorer.usgs.gov/). Michael Eineder is thanked for fruitful discussions and encouragements regarding the work on this paper. Martin Zimmer, Christian Kujawa, Elske de Zeeuw-van Dalfsen, Nicole Richter, Mehdi Nikkho, Jacky Salzer, Giancarlo Tamburello, and Ayleen Gaete for joint field work activities and discussions. The editors and anonymous reviewers are thanked for their comments that helped improving the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Estimation of PWV Contents in the Volcanic Cloud**

Water vapor contents in the volcanic plume were estimated by scaling the SO2 column density profiles measured by our scanning DOAS with the molar H2O/SO2 ratio derived from gas compositional measurements of the Multi-GAS instrument. For this purpose we used the 'raw' differential SO2 SCD output of the NOVAC-software [91], which in our case provides SO2 path length concentrations in units of ppm·m, that is SO2 concentration distributed along the effective light paths captured at each observing angle of the scans. Several conversions of the 'raw' data were necessary, in order to estimate the precipitable water vapor content of the volcanic cloud in a vertical atmospheric column right above each point of the plume cross-sections.

The light paths associated with the conical viewing geometry of the scanning DOAS instrument are generally inclined with respect to the vertical, entailing amplification of the absorption signal in relation to the signal that would be obtained along the more direct zenith light path due to extension of the distance that the incoming light has to travel through the absorbing air mass. Extension of the effective light path along inclined viewing directions is typically taken into account by means of the so-called air mass factor (AMF), which relates the measured SO2 SCDs to their respective vertical column densities (VCDs). The SO2 VCDs thus simply correspond to the ratio of SO2 SCDs and air mass factor as expressed in Equation (A1) (e.g., [92]).

$$SO\_2\ VCD\ (ppm\ \cdot m) = \frac{SO\_2\ SCD}{AMF\_G} \tag{A1}$$

Here, we used a geometrical approximation of the air mass factor *AMFG*, that can be derived for each viewing direction of a scan by means of Equation (A2) [91]

$$AMF\_G = \sqrt{\frac{\left(\cos(\beta)\cos(\delta) + \sin(\beta)\cos(\theta)\sin(\delta)\right)^2 + \left(\sin(\beta)\sin(\theta)\right)^2}{\left(\cos(\beta)\sin(\delta) + \sin(\beta)\cos(\theta)\cos(\delta)\right)^2} + 1} + 1,\tag{A2}$$

where *θ* is the scan angle (or zenith angle, which corresponds to the angle between observing elevation and zenith position), *δ* = 0 is the tilt of the scanner of our DOAS station at Láscar, and *β* = 60 is half the opening angle of the conical scanning surface. More sophisticated radiative transfer corrections were not performed additionally. The resulting SO2 VCDs are given as path length concentrations in units of ppm·m, and were used to approximate the corresponding PWV contents in units of mm water column. This required several steps of conversion, which were combined in Equation (A3).

$$PWV\ (mm) = 10\ mm \times \frac{SO\_2\ VCD(ppm\text{-}m) \times 2.5 \times 10^{15} \text{SO}\_2\ molecules\cdot\text{cm}^{-2} \times \text{H}\_2\text{O}/\text{SO}\_2}{3.34 \times 10^{22} \text{H}\_2\text{O}\text{ molecules\cdot cm}^{-2}}\tag{A3}$$

Single steps of this conversion equation were derived from the following relationships. A path length concentration of 1 *ppm*·*m SO*<sup>2</sup> roughly corresponds to a column density of 2.5 × 1015*SO*<sup>2</sup> *molecules*·*cm*−<sup>2</sup> at standard conditions for temperature and pressure, that is 298.15 K and 1 atmosphere air pressure (e.g., [91]). The corresponding number of vertically integrated water vapor molecules per square centimeter thus is readily obtained by scaling the SO2 VCDs in units of molecules·cm−<sup>2</sup> with the H2O/SO2 molar ratio determined from the Multi-GAS measurements. H2O VCDs in units of molecules·cm−<sup>2</sup> further can be converted to the total amount of water vapor present in a vertical atmospheric column, also known as column integrated water vapor (IWV), which is commonly stated as the vertically integrated mass of water vapor

per unit area (e.g., kg·m−2). At conditions of 1 atmosphere air pressure and 273.15 K one can find that an IWV of 1 *g H*2*O cm*−<sup>2</sup> = 10 *kg H*2*O m*−<sup>2</sup> roughly corresponds to a column density of 6.02214 × 1023 *molecules*·*mol*<sup>−</sup>1/18.015 *<sup>g</sup>*·*mol*−<sup>1</sup> = 3.34 × 1022 *<sup>H</sup>*2*O molecules*·*cm*<sup>−</sup>2, and respectively to 10 *mm* of an equivalent column of liquid water (A4) (pp.3 and 4 of the Appendix in McClatchey et al. [93]; or p.3 in Wagner et al. [94]).

$$13.34 \times 10^{22} H\_2O \text{ molecules} \cdot \text{cm}^{-2} = 1 \text{ g} \cdot \text{cm}^{-2} = 10 \text{ kg} \cdot \text{m}^{-2} \approx 10 \text{ mm water column} \tag{A4}$$

Conversion of IWV in units of kg·m−<sup>2</sup> to PWV in units of mm water column is obtained using Equation (A5) (e.g., [61])

$$PWV(mm) = \frac{IWW}{\rho\_{LW}},\tag{A5}$$

where the density of liquid water *<sup>ρ</sup>LW* typically roughly is equal to 1 *<sup>g</sup>*·*cm*<sup>−</sup>3, and thus a water vapor column density of 1 *kg*·*m*−<sup>2</sup> roughly corresponds to 1 *mm* height of an equivalent column of liquid water.

#### **Appendix B. Compensation of Downwind Evaporation**

LWCs in the plume above the crater rim were unknown, thus it was not possible to determine the actual evaporation rate from a predefined amount of liquid water. Considering increasing water vapor concentrations, respectively increasing H2O/SO2 ratios in the downwind portion of the volcanic plume, we thus calculated the drying power of the ambient air, respectively potential evaporation rates at plume height, based on climatic variables obtained from GDAS1 soundings, and terrain surface roughness. Potential evaporation rates were then used to scale our fixed H2O/SO2 ratio in order to compensate for downwind evaporation.

Potential evaporation rates correspond to the amount of liquid water that can be transformed to vapor through evaporation, if sufficient water is available (e.g., [95]). Assuming an unlimited availability of water, evaporation can in such a situation continue until the air above the evaporating surface approaches saturation humidity, and therefore actual evaporation cannot exceed potential evaporation. In the present case, availability of liquid water does not seem to be a limiting factor inside the plume, since we made the observation that the plume typically still is partly condensed when it arrives at the scanning DOAS. We thus assume that scaling by the potential evaporation rate yields a feasible approximation of the increase in water vapor at the expense of liquid water in the downwind portion of the volcanic cloud. We further assume that downwind dilution of water vapor due to dispersion of the plume is largely covered by the variations of the SO2 column densities measured by the DOAS, since both water vapor and SO2 are equally affected by diffusion and turbulent mixing processes (e.g., [96]). The amount of evaporation occurring inside the plume is largely controlled by the humidity of the entrained ambient air and the degree of turbulent mixing with air parcels of the plume, which is mainly governed by wind speed and atmospheric stability. Evaporation is thus commonly more pronounced during periods with dry atmospheric conditions and high wind velocities, which cause the transport of the plume to be more turbulent, and it is less significant, when the atmosphere is more humid and less turbulent, due to small transport velocities.

Evaporation rates were calculated by an aerodynamic (respectively mass-transfer) method (e.g., [97]), using temperatures, dew point temperatures, and wind speeds obtained from GDAS1 soundings provided by the National Oceanic and Atmospheric Administration (NOAA) [79]. The aerodynamic approach can be expressed by a Dalton-type Equation (A6) [98,99],

$$E\_P\left(\mathfrak{c}m\cdot\mathfrak{sec}^{-1}\right) = f(\mathfrak{u})(\mathfrak{e}\_{\mathfrak{s}} - \mathfrak{e}),\tag{A6}$$

in which *EP* is the rate of potential evaporation (cm·sec−1), *<sup>f</sup>*(*u*) is a wind speed function (cm·sec−1·hPa<sup>−</sup>1), respectively the vapor transfer coefficient, which is based on a logarithmic vertical wind velocity profile, and the second term (*es* − *e*) is the humidity term corresponding to the vapor pressure deficit of air (hPa). The product of both is typically referred to as the drying power of air. The wind function of the Dalton-type equation can be written as (A7) (Equation (4) in [100])

$$f(u)\left(cm\cdot sec^{-1}\cdot hPa^{-1}\right) = \frac{\rho\_{AIR}\cdot \varepsilon k^2}{\rho\_{LW}\,^P}\,\frac{u\_1}{\ln\left(z\_1/z\_0\right)^2},\tag{A7}$$

where *<sup>ρ</sup>AIR* and *<sup>ρ</sup>LW* correspond to the density of moist air and liquid water (g·cm−3), respectively. *-* = 0.622 is the ratio of the molar mass of water vapor (18.015 *<sup>g</sup>*·*mol*<sup>−</sup>1) to the average molar mass of dry air (28.9647 *<sup>g</sup>*·*mol*<sup>−</sup>1), and *<sup>k</sup>* = 0.4 is the dimensionless von Kármán's constant describing the logarithmic velocity profile of a turbulent air flow near a rough boundary with a no-slip condition [101], i.e., it is assumed that wind speed is zero directly above the surface. *u*<sup>1</sup> refers to the wind speed (cm·sec−1) at measurement height *<sup>z</sup>*<sup>1</sup> (cm above ground surface), and *<sup>P</sup>* is barometric pressure (mbar or hPa) at measurement height *z*1. *z*<sup>0</sup> is the surface roughness height (cm above surface), which corresponds to the average height of obstacles in the trajectory of the wind [102]. The density of moist air *<sup>ρ</sup>AIR* (*g*·*cm*−*3*) was calculated using the ideal gas law (A8)

$$
\rho\_{AIR} \left( \mathcal{g} \cdot \text{cm}^{-3} \right) = \frac{0.1P}{R\_d (1 + 0.608qv) (T + 273.15)},\tag{A8}
$$

where *Rd <sup>J</sup>*·*kg*−1·*K*−<sup>1</sup> = 287.04 is the gas constant for dry air, *P* is barometric pressure (hPa), *T* air temperature (◦C), and *qv* specific humidity (g·g<sup>−</sup>1), which in turn was approximated by Equation (A9).

$$q\upsilon\left(\mathbf{g}\cdot\mathbf{g}^{-1}\right) = \epsilon\frac{\mathcal{e}}{P} \tag{A9}$$

Here, *-* = 0.622 again is the ratio of the molar mass of water vapor to the average molar mass of dry air, *e* is the partial pressure of water vapor in hPa, and *P* is barometric pressure in hPa.

Saturation vapor pressure *es* (hPa) was calculated using the August–Roche–Magnus formula (A10) [103] with coefficients determined by Sonntag [104] (Table 1 in [105]).

$$
\sigma\_s \, (hPa) = \, 6.112 \exp\frac{17.62T}{243.12 + T} \,\tag{A10}
$$

where *T* is temperature in Celsius degree. Water vapor partial pressure *e* (hPa) was calculated as a function of relative humidity and saturation vapor pressure using Equation (A11)

$$
\varepsilon \left( hPa \right) = \frac{RH}{100} \varepsilon\_{s\_{\prime}} \tag{A11}
$$

where *RH*/100 is the fractional relative humidity, and *es* saturation vapor pressure at the ground level of the GDAS1 sounding. Relative humidity (%) was obtained by comparison of actual water vapor partial pressure *e* and saturation vapor pressure *es*, using Equation (A12) with coefficients determined by Alduchov and Eskridge [105].

$$RH\left(\%\right) = 100 \frac{e}{e\_s} = 100 \frac{\exp\left(17.625 T\_D\right) / \left(243.04 + T\_D\right)}{\exp\left(17.625 T\right) / \left(243.04 + T\right)}\,\tag{A12}$$

where *T* and *TD* are air temperatures, and respectively dew point temperatures in Celsius degree. Density of liquid water *<sup>ρ</sup>LW* (kg·m<sup>−</sup>3) was calculated from temperatures *<sup>T</sup>* ( ◦C) at plume height using Equation (A13) that was empirically determined by Jones and Harris [106].

$$\begin{array}{l} \rho\_{LW} \left( \text{kg} \cdot \text{m}^{-3} \right) = 999.85308 + 6.32693 \times 10^{-2} \, T - 8.523829 \times 10^{-3} \, T^2 \\ \quad - 6.943248 \times 10^{-5} T^3 - 3.821216 \times 10^{-7} T^4 \end{array} \tag{A13}$$

The roughness height *z*<sup>0</sup> was approximated by Equation (A14), which was experimentally determined by Plate and Quraishi [107] in wind tunnel experiments, and which yields fair results in the absence of more precise information [108].

$$z\_0 \, (cm) = 0.15h\_\prime \,\tag{A14}$$

where *h* is the average height of roughness elements (cm). As in this work the volcano is the main obstacle, and since the evaporating surface, respectively the volcanic plume typically is located roughly at summit altitude of the volcano, the average height of roughness elements relevant for the turbulent mixing of the plume with the atmosphere, was assumed to be equal to the height of the volcano above the surrounding plateau (about 750 m).

#### **Appendix C. Estimation of PWV Contents in the Atmosphere**

Background PWV contents in the atmosphere above the volcano were approximated by means of an empirically determined regression Equation (A15), which is based on the experimental results of Garrison and Adler [109] and proved to be applicable over a broad range of climatic conditions (e.g., [110]). The calculation requires climatic variables including pressure, dew point temperature and temperature, which were obtained from vertical atmospheric GDAS1 profiles provided by NOAA.

$$PVV(mm) = \frac{4.1173 \ RHP}{1013.25(273.15 + T)} e\_s + 0.2,\tag{A15}$$

where *RH* is relative humidity (%), which was calculated using Equation (A12), *P* is barometric pressure in hPa, *T* is temperature in Celsius degree, and *es* saturation vapor pressure at the ground level of the GDAS1 sounding was calculated using Equation (A10). Note that we reduced the intercept of the equation from +2 to +0.2 mm water column, in order to adapt the formula to hyper-arid conditions that are characterized by average PWV contents of less than 1 mm.

#### **Appendix D. Wind Field during SAR Acquisitions**

Wind, in particular over arid regions, plays a major role for the transport of humidity in the atmosphere, which is also true for the humidity emitted from volcanoes. In order to understand humidity variations in the atmosphere above a volcano, it is thus necessary to know about the local wind field. The temporal variations of the wind field above Láscar volcano were examined using atmospheric soundings of the Global Data Assimilation System (GDAS1) with one-degree grid spacing provided by NOAA. The soundings were concatenated to form a time series of wind speed and direction encountered at Láscar volcano during the period of interest (Figure A1a,b). Wind speeds at plume altitude ranged from 0.1 to 18 m·sec−<sup>1</sup> during the period of interest, and generally decreased throughout the whole atmospheric profile during austral summer (December to February) (Figure A1a). High altitude winds over the region usually are prevailing westerly throughout the year (orange colors in the upper part of Figure A1b), and are interrupted by short periods, which are dominated by easterly wind anomalies during austral summer. Surface winds typically are oriented according to topography, and exhibit pronounced diurnal variations (alternating blue and orange colors in the lower part of Figure A1b). These diurnal variations are characterized by valley winds (anabatic upslope breezes depicted by orange colors in the lower part of Figure A1b) that typically occur during daytime due to insolation of mountain flanks, and increasing wind speeds towards the evening (light blue colors in the lower part of Figure A1a), whereas rather calm mountain winds (katabatic downslope breezes and over-hill flows) prevail during the night (depicted by blue colors in the lower part of Figure A1b and dark blue colors in Figure A1a). Diurnal variations in surface wind directions are more prominent during summer, frequently resulting in a characteristic diurnal pattern of plume transport directions. We observed that in the morning hours of a summer day the plume typically drifts towards west or

northwest, turning over south towards east before noon, where it stays more or less steady until sunset. Downslope winds frequently force the plume to follow the steep morphology of the volcano.

Time series of wind directions in the range of typical plume heights, namely several 100 m above and below the summit of Láscar volcano were extracted from the 500 and 550 mbar pressure levels of the GDAS1 soundings, respectively, in order to narrow down the range of possible plume transport directions at the times of SAR acquisitions (Figure A1c,d). Wind directions at the time of SAR acquisitions of track 111, which were conducted around sunrise at about 10:00 a.m. UTC (Figure 2b), were predominantly westerly during spring 2013 and late summer 2014, whereas easterly directions prevailed during early summer 2013. Wind directions were generally very similar over the broad altitude range of both pressure levels, indicating relatively stable wind conditions. Pronounced differences in wind directions above and below summit were mainly observed during early summer, when winds tended to be less stable with respect to the direction, and were characterized by low velocities.

**Figure A1.** (**a**,**b**) Time series of vertical atmospheric profiles depicting variations in (**a**) wind speed and (**b**) wind direction. Note the weak winds accompanied by strong variations of wind directions during austral summer (ranging from December 2013 to February 2014). (**c**,**d**) Time series of wind directions some 100 m (**c**) above and (**d**) below the summit of Láscar (500 and 550 mbar pressure levels of the GDAS1 soundings, respectively). Wind directions at the time of SAR observations are indicated by *red* (master scene) and *blue dots* (slave scene). Wind directions during austral summer (ranging from December 2013 to February 2014) were predominantly easterly at the time of SAR observations.

We further complemented our observations by means of an additional high resolution wind field analysis, which is based on 3-dimensional gridded hindcasts of the Weather Research and Forecasting model (WRF). This we did, in order to resolve how the wind field evolved from turbulent flow close to the ground surface to laminar flow at higher altitude during the times of SAR acquisitions, because the gas plume typically is being emitted at the interface between surface flows and higher altitude flows. To this end, the WRF simulation was run using a 900 m horizontal grid spacing and a vertical division of 51 (terrain-following) eta levels.

Surface wind fields obtained from the lowest eta level of the WRF simulation, reveal that katabatic (downslope) winds prevailed during most SAR acquisition times (Figure A2).

**Figure A2.** Time series of surface wind fields obtained from the lowest eta level of WRF, determined for acquisition times of each SAR image. Katabatic (downslope) mountain winds prevail at the time of SAR acquisitions, which were recorded during the early morning hours at about 10:04 a.m. (UTC), respectively 07:04 a.m. (CLST).

The wind fields of all acquisition times were averaged, in order to match the cumulative SAR delay estimates produced by WBDD. This was done for the first eta level, which represents the wind conditions closest to the surface (Figure A3a). Changing conditions with increasing altitude were captured by means of gradually increasing the number of eta levels, which were included into the averaged wind field. In other words eta levels 1–5, 1–10, and 1–15 were combined to yield averaged wind fields (Figure A3b–d). Average wind direction of the terrain following winds close to the surface was predominantly towards southeast aloft the plateau, and heading downslope, which is mainly oriented towards west, over the western flank of the plateau (Figure A3a). Gradual integration of more eta levels into the average wind field nicely illustrates the turning direction of the winds

with increasing altitude (Figure A3a–d), finally arriving at upper tropospheric westerly trade wind directions, which are typical for this region (Figure A3d).

**Figure A3.** Averaged wind fields combining wind fields of all SAR observation times. (**a**) Average surface winds from eta level 1 (**b**) Average winds up to 387 m above surface (**c**) Average winds up to 1402 m above surface (**d**) Average winds up to 3193 m above surface.
