Radar Path Delay Effects in Volcanic Gas Plumes: The Case of Láscar Volcano, Northern Chile

<jats:p>Modern volcano monitoring commonly involves Interferometric Synthetic Aperture Radar (InSAR) measurements to identify ground motions caused by volcanic activity. However, InSAR is largely affected by changes in atmospheric refractivity, in particular by changes which can be attributed to the distribution of water (H2O) vapor in the atmospheric column. Gas emissions from continuously degassing volcanoes contain abundant water vapor and thus produce variations in the atmospheric water vapor content above and downwind of the volcano, which are notably well captured by short-wavelength X-band SAR systems. These variations may in turn cause differential phase errors in volcano deformation estimates due to excess radar path delay effects within the volcanic gas plume. Inversely, if these radar path delay effects are better understood, they may be even used for monitoring degassing activity, by means of the precipitable water vapor (PWV) content in the plume at the time of SAR acquisitions, which may provide essential information on gas plume dispersion and the state of volcanic and hydrothermal activity. In this work we investigate the radar path delays that were generated by water vapor contained in the volcanic gas plume of the persistently degassing Láscar volcano, which is located in the dry Atacama Desert of Northern Chile. We estimate water vapor contents based on sulfur dioxide (SO2) emission measurements from a scanning UV spectrometer (Mini-DOAS) station installed at Láscar volcano, which were scaled by H2O/SO2 molar mixing ratios obtained during a multi-component Gas Analyzer System (Multi-GAS) survey on the crater rim of the volcano. To calculate the water vapor content in the downwind portion of the plume, where an increase of water vapor is expected, we further applied a correction involving estimation of potential evaporation rates of water droplets governed by turbulent mixing of the condensed volcanic plume with the dry atmosphere. Based on these estimates we obtain daily average PWV contents inside the volcanic gas plume of 0.2–2.5 mm equivalent water column, which translates to a slant wet delay (SWD) in DInSAR data of 1.6–20 mm. We used these estimates in combination with our high resolution TerraSAR-X DInSAR observations at Láscar volcano, in order to demonstrate the occurrence of repeated atmospheric delay patterns that were generated by volcanic gas emissions. We show that gas plume related refractivity changes are significant and detectable in DInSAR measurements. Implications are two-fold: X-band satellite radar observations also contain information on the degassing state of a volcano, while deformation signals need to be interpreted with care, which has relevance for volcano observations at Láscar and for other sites worldwide.</jats:p>

section_content: Volcano monitoring is effectively based on multi-parametric datasets, often characterizing surface deformation, degassing, seismicity and temperature changes in time and space. However, even though comparison of independent datasets is common, they are only rarely integrated into a common evaluation scheme to enhance the relative strengths of the individual applied methods. Seismic and deformation data for instance are increasingly being integrated to improve spatio-temporal localization of magma movement beneath a volcano, which enables to better constrain the geometries of magma pathways and reservoirs, and to derive more precise physical models of the associated processes [1, 2]. Gas emission and deformation data, however, have not yet been deeply considered as integrated information, though it has been shown that variations in degassing and deformation are often intimately coupled [3]. A combined analysis of time series measurements of volcanic degassing and deformation rates proved to be a powerful tool for the evaluation of volumetric changes in magma reservoirs due to the accumulation or discharge of volcanic gases (e. g., [4, 5] ). Dynamics of degassing are furthermore critical for understanding pressure fluctuations in a magmatic and hydrothermal system, which additionally may manifest as temporary deformation of the volcanic edifice [6] [7] [8]. However, these datasets are commonly considered separately, and the fate of gas emissions downwind of the emission source is typically not taken into account, which is why the contribution of gas emissions on the radar signal used for deformation measurements has not been comprehensively demonstrated. Differential interferometric synthetic aperture radar (DInSAR) measurements allow for the detection of mm-scale line of sight (LOS) displacements of the observed surface [9]. It has been well demonstrated that the accuracy of satellite radar ground deformation measurements is affected by changes in atmospheric refractivity, in particular by changes which can be attributed to the highly variable distribution of water vapor in the observed atmospheric column [10, 11]. 

Atmospheric contributions to DInSAR data often have similar magnitudes and wavelengths as the actual ground deformation signal, and thus they need to be removed from interferograms, when deformation measurements are the purpose of monitoring. There have been two main strategies for this, (i) time-space-based filtering, and (ii) modeling of atmospheric contribution [12]. In particular, the latter is highly successful, as interfering contributions from a moist atmosphere can coarsely be predicted and compensated by means of high-resolution numerical weather models (e. g., MM5, and WRF), which provide the information required to predict the atmospheric phase delay in the LOS of an interferometric measurement [13] [14] [15] [16]. 

Degassing volcanoes; however, produce their own atmospheric disturbances, which are typically not well captured by weather models [17, 18]. The large and variable amounts of water vapor in volcanic plumes may cause differential phase errors in interferometric measurements due to reduction of radar propagation velocity within the plume above and downwind of the volcano [19] [20] [21] [22]. 

The leeward sector of a degassing volcano is hence prone to be affected by pronounced differential phase signatures, which can be misinterpreted as a deforming ground surface, and may thus obscure the real ground deformation information of an interferogram. Rosen et al. [19] described the possible consequences that such propagation delays would have on DInSAR measurements using the examples of phase signatures that occurred on Kilauea volcano, Hawaii. Similarly, Wadge et al. [20] observed enhanced tropospheric delays in radar data of Soufrière Hills. Instead of using weather models for the correction of their DInSAR data, these authors exploited data from Global Positioning System (GPS) measurements on the leeward side of the volcano. Based on local SO 2 flux and gas compositional measurements, and assuming a homogeneous distribution of water vapor inside the volcanic gas plume, they however obtained an estimate of the gas plume related radar delay of only about 0. 05 mm, which is far below, and thus negligible with respect to InSAR resolution. Gas plume related disturbances are nevertheless also identified elsewhere in interferograms of degassing volcanoes (e. g., at Pico do Fogo as described in González et al. [21] ), but remained to be quantitatively demonstrated and validated for a selected showcase. 

To reach this aim, the precipitable water vapor (PWV) content in the plume at the time of SAR acquisitions has to be determined. PWV is the amount of water vapor (in kg•m -2 or mm) vertically integrated in an atmospheric column. Space borne radar interferometric delay measurements in principle can be used to map the spatial distribution of water vapor in the atmosphere [23] [24] [25], if the contributions from topographic phase and ground deformation are known, or when independent PWV estimates are used as a reference. Accordingly, those measurements also enable to map the spatial extent of a volcanic gas plume when the strength of the degassing source is known, and independent PWV estimates are integrated into analysis. 

Even though water vapor by far is the largest and most abundant component in the gas emissions of most volcanoes, volcanic water vapor emissions to date only rarely have been measured and quantified by means of remote sensing (and preferentially by means of indirect methods as e. g., in [26] ), due to the missing contrast relative to the high water vapor concentrations in the ambient atmosphere, which commonly preclude its quantification in the volcanic cloud. During the present decade, however, some progress has been made in this respect, and volcanologists have adopted and adjusted several technologies known from meteorological applications to remotely map and quantify volcanic water vapor emissions despite a variable background atmosphere. Fiorani et al. [27] for instance, successfully measured water vapor fluxes of Stromboli volcano in Italy using a ground-based LiDAR system, and Bryan et al. [28] introduced a combined millimeter-wave radar/radiometer system (also ground-based), which enables 3-dimensional mapping of the spatial distribution of ash and water vapor in eruptive clouds. These methods have in common that they use active sensors, which transmit a signal that is able to penetrate the volcanic cloud and thus allow measuring the reflection, refraction, or scattering of the signal occurring inside the cloud. Such methods thus are largely superior to methods deploying passive optical imaging sensors (such as UV-and IR-spectroscopical methods), which typically are strongly affected by the presence of liquid water and aerosols [29], and therefore fail to work in very dense clouds. Ground-based active sensors are, however, rather exotic, and are cost-and labor-intensive, and thus have not yet routinely been used for permanent gas emission monitoring purposes. 

We thus used a combination of more common sensors, which enable unsupervised continuous monitoring. Because direct measurements of volcanic water vapor emissions by means of optical remote sensing are challenging [30], we measured SO 2 column density profiles in conjunction with molar H 2 O/SO 2 ratios in this work, in order to estimate the apparent PWV contents in the gas plume above a volcano. These estimates were integrated as an a-priori knowledge into the analysis of DInSAR observations in order to investigate the propagation delay that was generated by volcanic gas emissions. We demonstrate the repeated contribution of the volcanic gas plume to the radar delay in the interferometric measurements and thus determine the plume-induced delay variations which can be falsely interpreted as a deformation signal.
section_name: Introduction

section_content: The study was conducted at Láscar volcano, one of the most active volcanoes of the central Andes [31, 32], which has been identified as the second largest emission source of volcanic gases in Northern Chile [33]. The volcano is located east of the Salar de Atacama basin (23. 37 • S, 67. 73 • W), on the western margin of the Puna plateau (Figure 1a ), which is one of the driest areas on Earth. Much of the plateau around the 5592 m high Láscar volcano has an altitude of more than 4 km above sea level. Background atmospheric PWV over the region is very low most of the year, and generally less than 1 mm total water column [34]. Cloudless skies, an extremely low atmospheric water content, and the high altitude lead to an exceptionally high atmospheric transparency [35], which in general is conducive for remote sensing applications and hence for the purpose of our study. 

The volcanic edifice of Láscar comprises two truncated intersecting composite cones, termed Western and Eastern edifice hereafter, and hosts 5 nested craters, which are aligned along an ENE-WSW trending lineament [31, 36, 37]. During the past two decades Láscar exhibited several periods of cyclic Remote Sens. 2018, 10, 1514 4 of 44 dome growth and collapses, which were punctuated by occasional explosive vulcanian to plinian eruptions [38]. Recent activity is characterized by gradual subsidence [39, 40] and persistently strong fumarolic degassing from the remnants of a dacitic dome [41] located in the westernmost crater of the Eastern edifice (Figure 1b, c ). The dome remnants cover an area of about 40,000 m 2 (about 230 m in diameter) in the central depression of the 800 m wide and 400 m deep crater. Numerous high temperature fumaroles (about 300 • C) are distributed on top of the dome and the surrounding inner crater walls [41] [42] [43]. Several low temperature fumaroles and diffuse vent sites are located in the eastern craters of the Eastern edifice and on the southern crater rim, but the high temperature gases ascending from the dome and the fumaroles in the active crater are the main emission source [41]. These hot gas emissions produce a persistent thermal anomaly in the middle of the crater, which has been subject of numerous remote sensing studies using infrared satellite imagery (e. g., [44] [45] [46] [47] [48] ). Gas emissions from Láscar volcano are very moist in contrast to the dry atmospheric background. Water vapor emissions from Láscar account for more than 90% of its total gas emissions, approaching rates of several metric kilotons per day [33]. The persistently degassing, and largely unvegetated stratovolcano (Figure 1d ) thus provides ideal conditions to examine radar propagation delays caused by refractivity variations in volcanic gas. 

Remote Sens. 2018, 10, x FOR PEER REVIEW 4 of 45 strong fumarolic degassing from the remnants of a dacitic dome [41] located in the westernmost crater of the Eastern edifice (Figure 1b, c ). The dome remnants cover an area of about 40,000 m 2 (about 230 m in diameter) in the central depression of the 800 m wide and 400 m deep crater. Numerous high temperature fumaroles (about 300 °C) are distributed on top of the dome and the surrounding inner crater walls [41] [42] [43]. Several low temperature fumaroles and diffuse vent sites are located in the eastern craters of the Eastern edifice and on the southern crater rim, but the high temperature gases ascending from the dome and the fumaroles in the active crater are the main emission source [41]. These hot gas emissions produce a persistent thermal anomaly in the middle of the crater, which has been subject of numerous remote sensing studies using infrared satellite imagery (e. g., [44] [45] [46] [47] [48] ). Gas emissions from Láscar volcano are very moist in contrast to the dry atmospheric background. Water vapor emissions from Láscar account for more than 90% of its total gas emissions, approaching rates of several metric kilotons per day [33]. The persistently degassing, and largely unvegetated stratovolcano (Figure 1d ) thus provides ideal conditions to examine radar propagation delays caused by refractivity variations in volcanic gas.
section_name: Study Area

section_content: In this section, we will introduce a method that enables us to isolate and map gas plume related phase delay contributions in interferometric radar measurements by means of correlating a-priori information on the temporal variations of PWV contents in the volcanic gas plume with the signal strength variations at each range azimuth position of a DInSAR time series. Application of the method further requires involving several different prior constraints encompassing the temporal variations of all other processes that potentially contribute to total DInSAR phase in order to derive an estimate of the gas plume related interferometric signal in the DInSAR time series. In the following we will specify which data were used for our case study spanning the period October 2013 to February 2014 (later often simply referred to as the considered period), present the techniques used to derive the necessary a-priori knowledge, and describe how the estimation technique works.
section_name: Data and Methods

section_content: SO 2 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 SO 2 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. 

SO 2 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 SO 2 [51], and O 3 [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].
section_name: SO 2 Column Density Retrieval

section_content: Water vapor contents in the Láscar gas plume were obtained by scaling the SO 2 SCDs with the molar H 2 O/SO 2 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 SO 2, H 2 S and H 2, and a non-dispersive closed-path IR spectrometer for measuring CO 2 and H 2 O 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 H 2 O/SO 2 ratio obtained for the measurement period, which was used to scale the SO 2 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 SO 2 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 SO 2 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.
section_name: Estimation of Water Vapor Contents in the Láscar Plume

section_content: 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).
section_name: SAR Data and InSAR Methods

section_content: 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

where δ 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 Π -1 T S (≈ 6. 5) is a dimensionless conversion factor that is approximated using surface temperature T S 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

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 (PWV DIFF = PWV master -PWV slave ), and used to determine the differential SWD (dSWD) that such differential PWV contents theoretically would produce in the gas plumes of each interferogram.
section_name: Determination of Gas Plume Related Phase Delays

section_content: 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.
section_name: DInSAR Preparation

section_content: 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-CWT), which is being used in a wide range of applications in the field of signal and image processing [64]. The DT-CWT 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-CWT 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-CWT 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. 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.
section_name: DInSAR Decomposition

section_content: 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 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.
section_name: Priors Used for DInSAR Decomposition

section_content: 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.
section_name: Priors Used for DInSAR Decomposition

section_name: Results

section_content: 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. SO 2 SCDs measured within the gas plume on the days with TerraSAR-X acquisitions ranged from 2. 5 × 10 16 to 2. 8 × 10 18 molecules•cm -2 and averaged at 7. 2 × 10 17 molecules•cm -2. SO 2 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 CO 2 /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 H 2 O/SO 2 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 H 2 O/SO 2 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 SO 2 SCDs at each scan angle of the plume cross sections using the fixed molar H 2 O/SO 2 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 H 2 O 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 H 2 O/SO 2 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 SO 2 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.
section_name: PWV Contents in the Volcanic Gas Plume of Láscar

section_content: 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.
section_name: Gas Plume Related Phase Delays

section_content: 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).
section_name: Isolation of the Gas Plume Related Signal in DInSARs

section_content: 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.
section_name: Residuals of the Refractivity Related Phase Delay and Ground Deformation

section_content: 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). 

here sgn denotes the sign function and

corresponds to the amplitude of the gas plume estimate obtained from WBDD and → 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% (R 2 = 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.
section_name: Validation of the Gas Plume Signal Estimate

section_name: Discussion

section_content: 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].
section_name: Why the Distinction between Periods of Volcano Deformation and Enhanced Degassing is Important

section_content: 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] [77] [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 H 2 O 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.
section_name: How Methodological Limitations Can be Turned into Benefit

section_content: 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.
section_name: Phase Delay Effects above the Crater and Downwind of Láscar

section_content: 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. SO 2 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 H 2 O/SO 2 molar ratio of 34:1 obtained from Multi-GAS, in order to retrieve the corresponding H 2 O-fluxes (Figure 7a ). 

Assuming a constant H 2 O/SO 2 ratio and a constantly strong condensation of the emitted water vapor, the H 2 O (SO 2 ) emission rates ranged from 5 to 100 kt•day -1 (150 to 2900 t•day -1 ) and averaged at 12 kt•day -1 (350 t•day -1 ) during the considered period. Based on these assumptions we estimate that the total errors of our upscaled H 2 O emission rates likely were around -15 to +45% and thus even more skewed towards underestimation than the error envelopes of corresponding SO 2 -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 H 2 O and respectively 550 tons of SO 2 per day) with respect to the rest of the period (10 kilotons of H 2 O and respectively 300 tons of SO 2 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 SO 2 emission rates confirmed this suspicion and appropriately showed sharp increases over the course of each of these days (Figure 7a ). 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.
section_name: Gas Emission Rates from Láscar Volcano

section_content: 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 (R 2 = 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 SO 2 (H 2 O) 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 SO 2 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 SO 2 -fluxes is the dissolution of emitted SO 2 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 SO 2 emissions is generally estimated being on the order of merely a few percent of total emitted SO 2 per hour. This in turn leads to relatively long atmospheric lifetimes of SO 2 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.
section_name: Estimation of Background Atmospheric PWV Contents

section_content: 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 SO 2 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 SO 2 and increased hydrothermal H 2 O contributions of meteoric origin. The assumption that molar H 2 O/SO 2 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 H 2 O/SO 2 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.
section_name: Limitations of the Proposed Method

section_content: 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.
section_name: Advantages of the Proposed Method

section_content: 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 H 2 O emission rates ranged from 5 to 100 kt•day -1 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. 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.
section_name: Summary and Conclusions

section_content: 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. 

wind velocity profile, and the second term (e se) 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] )

where ρ AIR and ρ 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 g•mol -1 ) to the average molar mass of dry air (28. 9647 g•mol -1 ), and k = 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 1 refers to the wind speed (cm•sec -1 ) at measurement height z 1 (cm above ground surface), and P is barometric pressure (mbar or hPa) at measurement height z 1. z 0 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 ρ AIR (g•cm -3 ) was calculated using the ideal gas law (A8)

where R d J•kg -1 •K -1 = 287. 04 is the gas constant for dry air, P is barometric pressure (hPa), T air temperature ( • C), and qv specific humidity (g•g -1 ), which in turn was approximated by Equation (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 e s (hPa) was calculated using the August-Roche-Magnus formula (A10) [103] with coefficients determined by Sonntag [104] where T and T D are air temperatures, and respectively dew point temperatures in Celsius degree. Density of liquid water ρ LW (kg•m -3 ) was calculated from temperatures T ( • C) at plume height using Equation (A13) that was empirically determined by Jones and Harris [106]. 

ρ LW kg•m -3 = 999. 85308 + 6. 32693 × 10 -2 T -8. 523829 × 10 -3 T 2 -6. 943248 × 10 -5 T 3 -3. 821216 × 10 -7 T 4 (A13)

The roughness height z 0 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]. 

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).
section_name: Acknowledgments:

section_content: 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.
section_name: Appendix C. Estimation of PWV Contents in the Atmosphere

section_content: 4. 1173 RH P 1013. 25(273. 15

where RH is relative humidity (%), which was calculated using Equation (A12), P is barometric pressure in hPa, T is temperature in Celsius degree, and e s 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.
section_name: PWV(mm) =

section_content: 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 -1 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. 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 ). 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 ). 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 ). 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 ).
section_name: Appendix D. Wind Field during SAR Acquisitions

section_content: Delay maps of single event APSs were generated to mitigate non-repeating atmospheric disturbances from the gas plume estimate, and were used to infer the meteorological conditions during each SAR observation (Figure A4 ). These APS estimates are presented in the form of simple phase delay maps, where the scale indicates lengthening or shortening of the radar delay in units of millimeters. Such phase delay maps are snapshots of the meteorological situation, reflecting the spatial distribution of water vapor fields at the times of SAR acquisitions. Atmospheric disturbances aloft volcanoes are generally more pronounced and show more complex flow patterns on the lee side of the volcanic edifice [18]. This anisotropic distribution of turbulent atmospheric patterns can be attributed to the presence of volcanic gas plumes on the one hand, but also to orographic effects that govern the transport of moist air over mountainous terrain [17, 20]. Orographic effects comprise diabatic heating of air masses over insolated mountain flanks and orographic lifting of air masses that are pushed by the wind, causing upslope advection of moist air that is forced to rise following the steep topography, and to cool adiabatically causing an increase of the relative humidity.
section_name: E.1. Single Event APS Estimates

section_name: Appendix E. The Decomposed APS: Non-Repeating and Repeating Atmospheric Phase Delays

section_content: Delay maps of single event APSs were generated to mitigate non-repeating atmospheric disturbances from the gas plume estimate, and were used to infer the meteorological conditions during each SAR observation (Figure A4 ). These APS estimates are presented in the form of simple phase delay maps, where the scale indicates lengthening or shortening of the radar delay in units of millimeters. Such phase delay maps are snapshots of the meteorological situation, reflecting the spatial distribution of water vapor fields at the times of SAR acquisitions. Atmospheric disturbances aloft volcanoes are generally more pronounced and show more complex flow patterns on the lee side of the volcanic edifice [18]. This anisotropic distribution of turbulent atmospheric patterns can be attributed to the presence of volcanic gas plumes on the one hand, but also to orographic effects that govern the transport of moist air over mountainous terrain [17, 20]. Orographic effects comprise diabatic heating of air masses over insolated mountain flanks and orographic lifting of air masses that are pushed by the wind, causing upslope advection of moist air that is forced to rise following the steep topography, and to cool adiabatically causing an increase of the relative humidity. APS delay patterns were compared to wind directions and atmospheric PWV contents estimated for the times of SAR acquisitions to assure that the modeled and measured wind directions and PWV estimates are consistent with the location and strength of associated phase delay patterns. Wind directions determined for the summit region of Láscar were prevailing westerly during the early morning hours of the "dry season", whereas wind directions of the "wet season" were predominantly easterly (see Figure A1c,d, Figures A2 and A3 ). This is reflected by the spatial distribution of cloud shaped signatures, which are confined to the plateau east of the volcano during the "dry season", while pronounced cloud shaped signatures are mainly confined to the western flank of the volcano during early summer, indicating that moist air has been transported at low-altitude by easterly winds towards the edge of the plateau during the "wet season" (period comprising SAR acquisitions of 23 December 2013 to 14 January 2014). Wind directions are thus consistent with the spatial distribution of observed phase delay patterns and PWV estimates agree with the strength of these patterns. Enhanced humidity variations encountered during SAR acquisitions of early summer can further be ascribed to the SAR acquisition strategy. Space based SARs typically repeat their observations at the same local time, causing SAR observations of early summer to be more affected by atmospheric disturbances, since they are recorded later with respect to sunrise, due to variations in the length of the day (Figure 2b ). At Láscar volcano this effect is additionally enforced by more humid conditions that generally occur during that period.
section_name: E.1. Single Event APS Estimates

section_content: The phase delay estimates obtained for air temperature, air pressure and relative humidity priors comprise refractivity related phase contributions, which repeatedly occurred in all DInSARs of the time series and thus were not captured by the single event APS estimates. Such phase contributions therefore may contain residues of the stratification, which have not yet been removed through the coarse atmospheric correction that was performed prior to WBDD analysis utilizing the phase delay simulations obtained from the WRF, as well as repeating orographic effects, which may occur in multiple interferograms due to similar weather conditions. 

Phase delay patterns in the air temperature dependent phase screen (Figure 6e ) are asymmetrically distributed with respect to topography of the volcanic edifice, which can be ascribed to the distinct exposure of mountain flanks to sunlight. All SAR observations used in this study were made around sunrise, thus the western flank commonly lies in the shadow of the volcanic edifice, resulting in cooler air masses with a higher refractivity aloft the western flank, which produces a lengthening of the propagation path delay, whereas the air masses above the flat plain southeast of the volcano, the summit region, and the eastern flanks are subjected to diabatic heating due to a more pronounced insolation, which causes refractivity in the overlying air mass to be smaller due to higher air temperature. 

Phase delay patterns of the pressure related phase screen (Figure 6f ) indicate a lengthening of the delay over several confined steep-sloped areas, which are particularly exposed to westerly winds. The phase delay patterns of the temperature and pressure phase screens (Figure 6e,f) thus have the opposite direction, if compared to the phase delay patterns of the gas plume estimates (Figure 6c, d ), resulting in a partial cancellation of the plume related phase delay (Figure 6b ). Cancellation of the opposing phase screens hence reflect the negative dependence of volcanic gas emissions on barometric pressure and ambient temperature [77, 78]. This is further supported by the spatial distribution of phase delay patterns in the relative humidity related phase screen (Figure 6g ), which indicate enhanced relative humidity over exactly the same areas, which show a decreased temperature in the temperature phase screen (Figure 6e ). which locally may have resulted in substantial changes in surface altitude that occurred previous to the period considered here. This in turn may have given rise to phase differences in our SAR interferograms, where the ground surface geometries measured by SAR observations and modeled digital elevation are different. The amplitude of such topographical artifacts is largely controlled by the measurement geometry of the DInSAR measurement, and generally increases proportionally to the length of the spatial baseline. As these phase differences can become exceedingly large in interferograms with a long spatial baseline, we instead used very small baseline interferograms for our DInSAR decomposition analysis. 

The patterns in the phase delay estimate that we obtained for the spatial baseline prior are indeed intimately linked to topographic features (Figure 6h ), and thus very likely reflect the development of the land surface during the period following the SRT Mission. The most prominent feature of this phase delay estimate is a pronounced shortening of the phase delay in the summit region of Láscar and on the SW flank of Aguas Calientes, which may be attributed to the deposition of erupted material. Furthermore, a lengthening of the phase delay occurs along morphological depressions, which follow the steep southeastern flank of Láscar and the base of Aguas Calientes, and therefore may be a result of removal or compaction of sedimentary deposits. 

Table A1. Absolute values of predicted dSWDs, obtained from GDAS1 soundings, along with estimates of the mean measured APS amplitudes (= scale bar fractions of DInSARs) and RMSDs of phase delay amplitudes in DInSARs. of the volcanic gas plume is expected to increase the refractivity with respect to its surroundings. Negatively correlated phase delay patterns (depicted by cyan to blue colors in Figure A6a,b), which are distributed over the western flanks of Láscar and Aguas Calientes volcanoes, in contrast are incoherent in terms of the expected sign of the phase delay (Equation (1) in Smith and Weintraub, [66] ) and thus cannot be attributed to the gas plume. They rather indicate the contribution of another process (or other processes), which was (were) not mitigated from the gas plume estimate. Subset 02 was additionally processed without the SAR acquisition of 3 January 2014, in order to display the resulting reduction of the phase delay amplitude in the repeating plume related signal (Figure A6c ). PWV contents inside the plume and their respective differences were larger for SAR observations of subset 01, than they were for observations of subset 02 (Table 2 ), suggesting a stronger influence on the gas plume related phase delay signal of subset 01 (Figure A6a ), if compared to subset 02 (Figure A6b, c ). The amplitudes of the estimated phase delays displayed in Figure A6a -c are thus consistent with their respective input dSWD values (Table 3, column 4), and are also in accordance with humidity contrast between plume and atmosphere (Table 2 ). 

Remote Sens. 2018, 10, x FOR PEER REVIEW 39 of 45 detailed description of the wind field is given). Furthermore, this delay lengthening is coherent regarding the expected effect (Equation (1) in Smith and Weintraub, [66] ), because the enhanced water vapor content of the volcanic gas plume is expected to increase the refractivity with respect to its surroundings. Negatively correlated phase delay patterns (depicted by cyan to blue colors in Figure A7a,b), which are distributed over the western flanks of Láscar and Aguas Calientes volcanoes, in contrast are incoherent in terms of the expected sign of the phase delay (Equation (1) in Smith and Weintraub, [66] ) and thus cannot be attributed to the gas plume. They rather indicate the contribution of another process (or other processes), which was (were) not mitigated from the gas plume estimate. Subset 02 was additionally processed without the SAR acquisition of 3 January 2014, in order to display the resulting reduction of the phase delay amplitude in the repeating plume related signal (Figure A7c ). PWV contents inside the plume and their respective differences were larger for SAR observations of subset 01, than they were for observations of subset 02 (Table 2 ), suggesting a stronger influence on the gas plume related phase delay signal of subset 01 (Figure A7a ), if compared to subset 02 (Figure A7b, c ). The amplitudes of the estimated phase delays displayed in Figure A7a -c are thus consistent with their respective input dSWD values (Table 3, column 4), and are also in accordance with humidity contrast between plume and atmosphere (Table 2 ).
section_name: E.2. Repeating Atmospheric Phase Delays

section_name: Conflicts of Interest:

section_content: The authors declare no conflict of interest.
section_name: Conflicts of Interest:

section_content: Water vapor contents in the volcanic plume were estimated by scaling the SO 2 column density profiles measured by our scanning DOAS with the molar H 2 O/SO 2 ratio derived from gas compositional measurements of the Multi-GAS instrument. For this purpose we used the 'raw' differential SO 2 SCD output of the NOVAC-software [91], which in our case provides SO 2 path length concentrations in units of ppm•m, that is SO 2 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 SO 2 SCDs to their respective vertical column densities (VCDs). The SO 2 VCDs thus simply correspond to the ratio of SO 2 SCDs and air mass factor as expressed in Equation (A1) (e. g., [92] ). 

Here, we used a geometrical approximation of the air mass factor AMF G, that can be derived for each viewing direction of a scan by means of Equation (A2) [91] AMF G = (cos(β) cos(δ) + sin(β) cos(θ) sin(δ)

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 SO 2 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). 

Single steps of this conversion equation were derived from the following relationships. A path length concentration of 1 ppm•m SO 2 roughly corresponds to a column density of 2. 5 × 10 15 SO 2 molecules•cm -2 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 SO 2 VCDs in units of molecules•cm -2 with the H 2 O/SO 2 molar ratio determined from the Multi-GAS measurements. H 2 O VCDs in units of molecules•cm -2 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

Conversion of IWV in units of kg•m -2 to PWV in units of mm water column is obtained using Equation (A5) (e. g., [61] )

where the density of liquid water ρ LW typically roughly is equal to 1 g•cm -3, and thus a water vapor column density of 1 kg•m -2 roughly corresponds to 1 mm height of an equivalent column of liquid water.
section_name: Appendix A. Estimation of PWV Contents in the Volcanic Cloud

section_content: 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 H 2 O/SO 2 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 H 2 O/SO 2 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 SO 2 column densities measured by the DOAS, since both water vapor and SO 2 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],

in which E P is the rate of potential evaporation (cm•sec -1 ), f (u) is a wind speed function (cm•sec -1 •hPa -1 ), respectively the vapor transfer coefficient, which is based on a logarithmic vertical
section_name: Appendix B. Compensation of Downwind Evaporation

section_content: The SRTM-1 DEM that we used as a reference surface for our DInSAR observations is based on data which have been recorded in
section_name: Appendix F. DEM Error Related Phase Delays

section_content: The atmospheric phase contributions in our coarsely corrected DInSAR maps (Figure 5b ) are mainly caused by the presence of water vapor. Amplitudes of the atmospheric contribution to total DInSAR phase thus to a certain extent may be predicted using the modeled atmospheric PWV contents that were obtained from GDAS1 soundings for the summit region of Láscar volcano (Figure 7b and Table 2, columns 6 and 11). 

In order to assess how well these modeled PWV values are suited to reproduce the APS amplitudes we therefore utilized the calculations described in Section 3. 4 to determine the differential slant wet delays (dSWDs) that the PWV contents would produce in respective DInSARs, and compared their absolute values to the measured amplitudes of associated phase delays in each of the corresponding DInSAR maps (Figure 5b ). 

As a first step, we compared these predicted dSWDs (Table A1, column 4) to visually determined estimates of the mean APS amplitude in each of the DInSARs (Table A1, column 5). Using the scale bar ranges of the DInSAR maps in Figure 5 as reference, the mean APS amplitudes in five of the DInSARs (DInSARs 1-3, 5 and 6) were estimated to roughly amount to half of the scales, because a majority of the APS amplitudes does not exceed half of the associated scale, while they rather span 2/3rds of the scale bar in DInSAR 4, and only about 1/4 in DInSAR 7. Comparing the absolute values of the predicted dSWDs with our estimates of the mean APS amplitudes measured in corresponding DInSARs reveals that the latter are surprisingly well represented by the modeled dSWDs (R 2 = 0. 96; Figure A5a ), which at average are offset from measured APS amplitudes by merely 2. 1 mm. 

In order to obtain more representative values for the measured amplitudes of the atmospheric delays in the DInSAR maps (Figure 5b ) we further determined the root-mean-square-deviation (RMSD) of the phase delay amplitudes in each DInSAR map (Table A1, column 6) and compared it to the absolute values of the dSWDs, which the GDAS1 soundings predicted for the acquisition times of corresponding DInSARs. The comparison revealed that the absolute values of the predicted dSWDs are in good agreement with the RMSDs of DInSAR amplitudes (R 2 = 0. 70; Figure A5b ).
section_name: Appendix G. Estimation of APS Amplitude using Modeled Atmospheric PWV Contents Obtained from GDAS1 Soundings

section_content: Absolute values of predicted dSWDs, obtained from GDAS1 soundings, along with estimates of the mean measured APS amplitudes (= scale bar fractions of DInSARs) and RMSDs of phase delay amplitudes in DInSARs.
section_name: Table A1

section_content: Interferograms of DInSAR time series subsets 01 and 02 were decomposed, in order to separately examine the associated gas plume related phase delays. Estimates of gas plume related phase delays were computed for each of the two disjoint time series subsets, i. e., using a limited number of three (subset 01), and respectively four temporally interconnected interferograms (subset 02). Presented are the results obtained from the WBDD run, that was conducted omitting surface temperature and pressure priors, i. e., the phase delays associated to refractivity changes caused by atmospheric pressure and temperature variations were not prevented from leaking into the phase delay estimates obtained for the gas plume. Gas plume estimates of both subsets show a fan-shaped interferometric pattern indicating lengthening of the radar path over the south-eastern flank of Láscar volcano (depicted by yellow to orange colors in Figure A7a, b ), which is in good agreement with the common direction of plume transport (see supplementary Figures A1c, d
section_name: Appendix H. Gas Plume Related Phase Delays of DInSAR Time Series Subsets 01 and 02

section_content: Interferograms of DInSAR time series subsets 01 and 02 were decomposed, in order to separately examine the associated gas plume related phase delays. Estimates of gas plume related phase delays were computed for each of the two disjoint time series subsets, i. e., using a limited number of three (subset 01), and respectively four temporally interconnected interferograms (subset 02). Presented are the results obtained from the WBDD run, that was conducted omitting surface temperature and pressure priors, i. e., the phase delays associated to refractivity changes caused by atmospheric pressure and temperature variations were not prevented from leaking into the phase delay estimates obtained for the gas plume. Gas plume estimates of both subsets show a fan-shaped interferometric pattern indicating lengthening of the radar path over the south-eastern flank of Láscar volcano (depicted by yellow to orange colors in Figure A6a, b ), which is in good agreement with the common direction of plume transport (see supplementary Figure A1c,d, Figures A2 and A3 in Appendix D, where a detailed description of the wind field is given). Furthermore, this delay lengthening is coherent regarding the expected effect (Equation (1) in Smith and Weintraub, [66] ), because the enhanced water vapor content
section_name: Appendix H. Gas Plume Related Phase Delays of DInSAR Time Series Subsets 01 and 02