Investigation of Potential Volcanic Risk from Mt. Baekdu by DInSAR Time Series Analysis and Atmospheric Correction

Mt. Baekdu is a volcano near the North Korea-Chinese border that experienced a few destructive eruptions over the course of its history, including the well-known 1702 A.D eruption. However, signals of unrest, including seismic activity, gas emission and intense geothermal activity, have been occurring with increasing frequency over the last few years. Due to its close vicinity to a densely populated area and the high magnitude of historical volcanic eruptions, its potential for destructive volcanic activity has drawn wide public attention. However, direct field surveying in the area is limited due to logistic challenges. In order to compensate for the limited coverage of ground observations, comprehensive measurements using remote sensing techniques are required. Among these techniques, Differential Interferometric SAR (DInSAR) analysis is the most effective method for monitoring surface deformation and is employed in this study. Through advanced atmospheric error correction and time series analysis, the accuracy of the detected displacements was improved. As a result, clear uplift up to 20 mm/year was identified around Mt. Baekdu and was further used to estimate the possible deformation source, which is considered as a consequence of magma and fault interaction. Since the method for tracing deformation was proved feasible, continuous DInSAR monitoring employing upcoming SAR missions and advanced error regulation algorithms will be of great value in monitoring comprehensive surface deformation over Mt. Baekdu and in general world-wide active volcanoes.


Introduction
Monitoring of the circum-Pacific tectonic activity and volcanism is paramount for understanding mechanisms and coupling of processes and, more importantly, for assessing and potentially mitigating risks imposed to human settling. For these reasons, governmental bodies have been putting resources in establishing in situ ground-based as well as orbital monitoring platforms. When compared to countries facing a long history of exposure to volcanic risk in Asia, such as Japan Observation of abnormally high thermal activity in hot springs.

5.
Record of surface deflation indicating new magma activity at the western and northern slopes beginning in 2009.
Consequently, the ground observations in 2010 identified the incensement of the potential volcanic activity over Mt. Baekdu. Recently, Choi et al. [5] conducted a gravity anomaly analysis using the Earth Gravitation Model (EGM) dataset covering Mt. Baekdu, and attempted to model the magma source. Due to Mt Baekdu's close vicinity to a densely populated area, the record of potential volcanic activity has drawn considerable attention. It is, therefore, of importance to develop monitoring techniques allowing to conduct long-term observations of Mt. Baekdu. After the initial research by Zebker and Goldstein [6], volcano monitoring with Synthetic Aperture Radar (SAR) data has been conducted by a number of researchers in subsequent years [7][8][9][10][11][12]. As DInSAR techniques have commonly been considered as one of the most favourable and effective remote-sensing solutions for monitoring of volcanic activities [13,14], we here conducted a DInSAR survey over Mt. Baekdu to assess its potential activity.
To do so, a number of technical challenges needed to be solved, including: 1.
Corrections due to dense vegetation canopy at the summit of Mt. Baekdu which creates low-phase coherence for the DInSAR analysis. 2.
Atmospheric corrections due to significant tropospheric errors over Mt. Baekdu.

3.
Sparse data in particular regarding topographic and atmospheric field and validation data.
Considering these challenges, a high-accuracy deformation detection scheme is required in order to evaluate the risk of volcanic eruptions. To address these limitations and confirm the feasibility of a DInSAR analysis for Mt. Baekdu, we tested two different approaches in this study.
Firstly, we employed a two-pass DInSAR analysis with L-band ALOS PALSAR, which was expected to produce robust observations in challenging environments due to its operation at long wavelengths. Analysis was conducted using atmospheric error correction by applying an external error-source map derived from a weather forecasting model. Secondly, SAR time series analysis techniques were applied using stable DInSAR analysis methods. These suppress atmospheric and topographic errors in surface deformation estimates by additionally employing C-band ENVISAT ASAR images.
In this study, the main focus was therefore put on establishing a reliable deformation analysis using a DInSAR methodology by minimizing the noise, and accounting for atmospheric errors.
including the localized water vapor component introduced by tropospheric turbulences and the vertical stratification. At last, as stated by Samsonov [27], the degradation of base Digital Terrain Model (DTM) quality may also contribute significant errors in DInSAR products. Hanssen [28] estimated a 1 cm deformation error by 1 m height difference in base DTM when performing DInSAR analysis. Hence the topographic error should be considered as well. . The faults lines and the estimated magma chambers from the gravity model in [5] and epicenters of earthquakes [2] are also represented. The Sentinel-2 space-borne image taken in October, 2010 is shown in (c). The snow cover, dense vegetation canopy and steep slope around the Cheonji Candela Lake can be observed.
Based on the challenges introduced above, we estimate the contribution of atmospheric and topographic errors. Initial investigation using Moderate Resolution Imaging Spectroradiometer (MODIS) data was first performed to understand the characteristics of the atmospheric error component. In Figure 2, the water vapor distribution on the 14th of April and 5th of May in 2009 is shown in the MODIS Total Perceptible Water (TPW) MOD 5 [29] products. The depletion of water vapor around the summit of Mt. Baekdu are obvious in Figure 2a and identified less clearly in Figure  2b. The orographic effect caused by precipitation on the flanks of high-altitude mountains or plateaus and the local dryness along the other side of the mountain are clearly identified at Mt. Baekdu. The strength of this effect partly depends on some climate factors such as wind direction, so that the amount of atmospheric error of DInSAR by orthographic effect varies temporally and is not easily predictable without external water vapor observation. Meanwhile, the refractivity in wave propagation by the stratified pressure in high altitude topography produces another atmospheric error in DInSAR analysis. Thus, these atmospheric related error factors over Mt. Baekdu should be taken into account. As for the topographic error, however, the ground and aerial surveying campaigns for geodetic control of topographic products are unavailable over the area of Mt. Baekdu. As a result, the investigation and interpretation of the deformation signal including the analysis of atmospheric error are the main technical topic of this study. . The faults lines and the estimated magma chambers from the gravity model in [5] and epicenters of earthquakes [2] are also represented. The Sentinel-2 space-borne image taken in October, 2010 is shown in (c). The snow cover, dense vegetation canopy and steep slope around the Cheonji Candela Lake can be observed.
Based on the challenges introduced above, we estimate the contribution of atmospheric and topographic errors. Initial investigation using Moderate Resolution Imaging Spectroradiometer (MODIS) data was first performed to understand the characteristics of the atmospheric error component. In Figure 2, the water vapor distribution on the 14th of April and 5th of May in 2009 is shown in the MODIS Total Perceptible Water (TPW) MOD 5 [29] products. The depletion of water vapor around the summit of Mt. Baekdu are obvious in Figure 2a and identified less clearly in Figure 2b. The orographic effect caused by precipitation on the flanks of high-altitude mountains or plateaus and the local dryness along the other side of the mountain are clearly identified at Mt. Baekdu. The strength of this effect partly depends on some climate factors such as wind direction, so that the amount of atmospheric error of DInSAR by orthographic effect varies temporally and is not easily predictable without external water vapor observation. Meanwhile, the refractivity in wave propagation by the stratified pressure in high altitude topography produces another atmospheric error in DInSAR analysis. Thus, these atmospheric related error factors over Mt. Baekdu should be taken into account. As for the topographic error, however, the ground and aerial surveying campaigns for geodetic control of topographic products are unavailable over the area of Mt. Baekdu. As a result, the investigation

Processing Strategy and Data Sets
To address atmospheric error issues, the water vapor correction method developed by Li et al. [30][31][32] was employed in this study. Since Medium Resolution Imaging Spectrometer (MERIS) sensor which was the source of water vapor information shares the same satellite platform with the ENVISAT ASAR system, there is no temporal gap exists between each SAR image acquisition and the water vapor data, allowing for reliable error estimates could be achieved [31,32].
In addition, methods based on InSAR time series calculation based on stacked interferograms allow to estimate and reduce the errors, such as Permanent Scatterers (PS) by Ferretti et al. [33,34] and SBAS by Berardino et al. [35]. However, although these methods are feasible for solving atmospheric error, the processing efficiency and required number of SAR image sequences become limitations in employing these approaches. Considering both the estimated small deformation and the low phase coherence in the target area, two-pass DInSAR processing using error correction based on a high-resolution weather forecasting model and improved interferogram time series analysis, were applied and cross-validated.

Two Pass DInSAR with Atmospheric Error Compensation
With the characteristic of long wavelength, the L-band SAR has significant merits for producing coherence when observing high altitude volcanoes which are covered with snow and vegetation. Therefore, in Mt. Baekdu target area, it is preferable to employ ALOS PALSAR which is currently the only available L-band spacebone SAR at the time of this paper preparation. However there are only16 SAR images covering three years time interval over Mt. Baekdu which is not suggested to establish time series in order to conduct stable PS and/or SBAS processing [36]. Thus, instead of time series analysis, we track the surface deformation using two pass ALOS PALSAR DInSAR analyses were conducted. Moreover, the terms involving electromagnetic wave propagation delay in troposphere over Mt. Baekdu were estimated and removed. The amount of tropospheric phase change was computed based on the formula using pressure estimation and corresponding height value and total water vapor. It is simply summarised as: where ZDDD is zenith dry delay difference and ZWDD is zenith wet delay difference, both in unit of , and is the incidence angle of SAR images. Basically the ZDDD can be estimated as the pressure and temperature which are correlated with the altitude. Thus, with a number of weather observations, the estimation and the correction of dry component can be achieved using the relationship listed below.

Processing Strategy and Data Sets
To address atmospheric error issues, the water vapor correction method developed by Li et al. [30][31][32] was employed in this study. Since Medium Resolution Imaging Spectrometer (MERIS) sensor which was the source of water vapor information shares the same satellite platform with the ENVISAT ASAR system, there is no temporal gap exists between each SAR image acquisition and the water vapor data, allowing for reliable error estimates could be achieved [31,32].
In addition, methods based on InSAR time series calculation based on stacked interferograms allow to estimate and reduce the errors, such as Permanent Scatterers (PS) by Ferretti et al. [33,34] and SBAS by Berardino et al. [35]. However, although these methods are feasible for solving atmospheric error, the processing efficiency and required number of SAR image sequences become limitations in employing these approaches. Considering both the estimated small deformation and the low phase coherence in the target area, two-pass DInSAR processing using error correction based on a high-resolution weather forecasting model and improved interferogram time series analysis, were applied and cross-validated.

Two Pass DInSAR with Atmospheric Error Compensation
With the characteristic of long wavelength, the L-band SAR has significant merits for producing coherence when observing high altitude volcanoes which are covered with snow and vegetation. Therefore, in Mt. Baekdu target area, it is preferable to employ ALOS PALSAR which is currently the only available L-band spacebone SAR at the time of this paper preparation. However there are only16 SAR images covering three years time interval over Mt. Baekdu which is not suggested to establish time series in order to conduct stable PS and/or SBAS processing [36]. Thus, instead of time series analysis, we track the surface deformation using two pass ALOS PALSAR DInSAR analyses were conducted. Moreover, the terms involving electromagnetic wave propagation delay in troposphere over Mt. Baekdu were estimated and removed. The amount of tropospheric phase change was computed based on the formula using pressure estimation and corresponding height value and total water vapor. It is simply summarised as: where ZDDD is zenith dry delay difference and ZWDD is zenith wet delay difference, both in unit of λ, and Θ is the incidence angle of SAR images. Basically the ZDDD can be estimated as the pressure and temperature which are correlated with the altitude. Thus, with a number of weather observations, the estimation and the correction of dry component can be achieved using the relationship listed below.
where ZDD m is zenith dry delay for the master image (see [37] for detailed model), ZDD s is zenith dry delay for the slave image, Ps (in Hpa) is the surface pressure, lat is the latitude (in degrees) of target area and H (in km) is the height. Compared with the dry delay correction, the effect of wet delay induced by water vapor is relatively significant and difficult to estimate. As a result, many different approaches have been proposed by which to estimate and rectify the magnitude of the wet delay. The model shown in [32] was used in this study: where ZWD m is zenith wet delay in master, ZWD s is zenith wet delay in slave, P w is the density of water, IWV is the total amount of water vapor over an observation point and c is a conversion factor for the total perceptible water vapor (TPWV) to ZWD. According to [38], c is very close to 6.67 although there is some dependency on temperature.
Since it is not feasible to achieve global water vapor distribution with dense enough GPS delay observations or radiometer measurements, Li et al. [31,32] used water vapor data derived from MODIS and MERIS as another solution. However, the cloud coverage issue and the temporal gap between each SAR image acquisition and the water vapor estimation limits the feasibility of the method. Although high accuracy water vapor data from MERIS with sufficient spatial resolution is available occasionally, the temporal gap between MERIS and ALOS PALSAR image acquisitions cannot be avoided. Especially, considering the high mountainous areas are seldom free from cloud coverage over the time acquisition of MERIS and MODIS, it is required to use another source of water vapor data.
In this study, we use high resolution water vapor maps derived from weather forecasting models to correct for two pass DInSAR in the error correction procedure, based on previous research conducted by Wedge et al. [39], Foster et al. [40] and Nico et al. [41]. [33] and SBAS [35] are two methods widely used for performing time series analysis. Nevertheless, a number of limitations were encountered when applying the two methods over the target area. For PS-InSAR processing, the temporal and spatial baseline conditions between SAR images in this study were inadequate for extracting persistent scatterers. This is because a minimum of 25 pairs is necessary to build reliable PS results [36]. In addition, because vegetation cover exists in the high-altitude area around the summit, it was not expected that a sufficient number of persistent scatterers could be extracted. The SBAS processing employs all the appropriate DInSAR pair with small baselines [35]. Once the stacked differential interferogram is constructed, a linear system consisting of a small baseline combination matrix, phase values, and mean phase velocities are calculated. Singular Value Decomposition (SVD) together with LP (low pass) and HP (high pass) filters are used to extract the spatial and temporal components of the noise. The atmospheric, orbital and base DEM artifacts are estimated and corrected, and then the deformation velocities can be derived. Although in common cases the SBAS approach is capable of achieving high density deformation rate, some interferograms with low phase coherence may produce unreliable results as a consequence of permanent snow and cover vegetation on Mt. Baekdu.

PS-InSAR
To resolve all these issues, a complementary approach developed by Hooper [42] for combining PS-InSAR and SBAS was employed in this study. Although conventional PS-InSAR results showed satisfactory performance for monitoring the target area where artificial scatterers were widely populated, as shown by Colesanti et al. [36], it is difficult to define the sufficient number of PS over the vegetated mountain area where the stable scatterer such as artificial structures are absent. Hooper et al. [43] tackled this problem by introducing a new method for PS definition by investigating spatial and temporal correlation of phase, and extracted reliable results over natural landscape terrain. The algorithm was later updated using spatial and temporal correlation of interferogram phase to find pixels with low-phase variance in all terrains [44] and produced reliable deformation over Volcano Alcedo (Galapagos) where the vegetation extended even up to the highest altitude. Subsequently, the algorithm was improved by incorporating the PS and SBAS results by applying the spatial and temporal correlation tracking method of phase for SBAS interferograms [42]. Consequently, this approach increased the spatial sampling resolution of deformation measurement, as well as the reliability of the phase cycle.
The measurement of PS pixel was conducted by measuring the variation of residual phase in the time series [43]: where N is the number of interferogram, Φ x,i is the phase change in pixel x, Φ x,i is the phase change in sample mean within window patch andΦ ex,i is the estimated phase change residual with window patch. If the value is smaller than a predefined threshold, it was considered as a stable scatterer. The StaMPS/MTI (Stanford Method for Persistent Scatterers/Multi-Temporal InSAR) software implementing InSAR persistent scatterer (PS) method and SBAS [42] was employed to determine the PS pixels over the target area with high challengeable conditions for PS tracking. Based on the identified time series PS pixels, the displacement occurring over the target areas could be identified. Successful monitoring in similar condition and approach were reported in Hooper et al. [44], Pinel et al. [45], Chang et al. [46] and Decriem et al. [47].
The StaMPS/MTI processing began with two pre-processing steps. The first was to focus the raw data, and the second was to form interferograms from single-look complex (SLC) images in which JPL's Repeat Orbit Interferometry PACkage (ROI-PAC) was applied for interferogram calculation [48]. Based on the coherence and amplitude data produced by ROI-PAC, PS candidates were selected and permanent PS pixels were further identified in the StaMPS routine. Then the Statistical-cost, Network-flow Algorithm for Phase Unwrapping (SNAPHU) package [49] was used to unwrap the PS phase. At the last stage, the StaMPS/MTI was implemented again to perform DEM error correction and atmospheric filtering. The ground displacement was also solved at this stage.

Data Sets
A number of 4 ALOS PALSAR images listed in Table 1 are used for two pass DInSAR processing. Additionally, in order to examine the reliability of the water vapor map produced by spaceborne data, the water vapor products from MODIS Total Perceptible Water (TPW) MOD 5 [29] was also acquired. Due to the consistent cloud coverage over the summit in Mt. Baekdu, it was judged visually that it was not feasible to find complete MODIS water vapor product covering target area. To overcome the limitation, Kriging interpolation capable of reconstructing natural 3D surface from largely voided spatial measurements was employed to spatially interpolate water vapor from MODIS TPW products. Then two water vapor maps before and after each SAR acquisition time are linearly temporally interpolated into the water vapor map for each SAR acquisition time. The resulting spatially and temporally interpolated water vapor maps for the first ALOS PALSAR DInSAR pair are shown in Figure 3. The water vapor map for the slave SAR images (Figure 3b) show proper water vapor depletion but the one for the master SAR image ( Figure 3a) shows a very different pattern from the expected orthographic water vapor distribution probably caused by an over-populated cloud cover. As the MODIS/MERIS water vapor products for another ALOS PALSAR pair possess even more cloud coverage and long temporal gaps, it was not possible to utilize the MODIS/MERIS water vapor product for the two pass DInSAR correction over Mt. Baekdu. Instead, the water vapor map produced through weather forecast modelling was applied.  ENVISAT ASAR imagery was chosen for the time series analysis owing to constant image acquisitions during the 2007-2010 period. A total of 19 ENVISAT ASAR C-band images in descending mode are available (see Table 2). It is sufficient to apply the proposed interferogram time series analysis in target area.  ENVISAT ASAR imagery was chosen for the time series analysis owing to constant image acquisitions during the 2007-2010 period. A total of 19 ENVISAT ASAR C-band images in descending mode are available (see Table 2). It is sufficient to apply the proposed interferogram time series analysis in target area.

WRF Model Processing
For the reasons stated above, the weather forecasting model was employed to create the water vapor and surface pressure products for atmospheric correction of DInSAR processing. The Fifth-Generation Penn State/NCAR Mesoscale Mode (MM5) (see [40,50]) used in precedent researches [41] successfully demonstrated the capability of the Weather Research and Forecasting (WRF) model for the error migration of InSAR. Hence the version 3.2.1 Advanced Research WRF (WRF-ARW) [51] was applied for the construction of weather model. The predicted domain consists of 60 × 60 grids with 2 km spatial resolution and 40 vertical levels covering ±3 h of corresponding SAR image acquisition time. For the initial and boundary condition of the WRF model, the National Center for Environmental Prediction Final (NCEP_FNL) Global Analysis data was used as a reference. The NCEP_FNL data have 1 • ×1 • grids with 26 vertical levels, and are produced every 6 h [52]. The microphysics was dealt with using the WRF Single-Moment 6-class (WSM 6) scheme [53]. The radiation was calculated using the Rapid Radiative Transfer Model (RRTM) scheme [54] for longwave radiation, and the MM5 scheme [55] for shortwave radiation. The surface layer was simulated using the MM5 similarity scheme [56][57][58]. The land surface was processed using the Noah and Surface Model [59]. The planetary boundary layer and urban surface were performed using Yonsei University (YSU) scheme [60] and Kain-Fritsch scheme [61], respectively.

Validation of the Atmospheric Error Model and Two Pass DInSAR Analyses Using ALOS PALSAR
Following the method and data described in Section 3.4, the WRF simulation was performed and the atmospheric models on the 4 ALOS PALSAR image acquisition dates were extracted. All quantities derived from WRF were subsequently converted to total water vapor following Equation (7): where G is the universal gravitation constant, q i is the water vapor mixture ratio in ith layer and p i is the pressure in ith layer. The computed result was used for the DInSAR wet delay correction. Once the two WRF maps for the master and slave images were constructed (Figure 4a,b), the ZWDD of DInSAR was extracted using Equation (4). It was clearly showed that the constructed water vapor from WRF successfully caught the water vapor void around the mountain summits due to orographic effects. Although the direct inter-comparison was not possible due to the temporal gap between MODIS water vapor and WRF simulation, the interpolated MODIS water vapor from the cloud free raw products and water vapor from WRF show similarities to some extents as shown in Figures 3b and 4b. Together with ZDD, based on the surface pressure in WRF and Equation (6), the total delay map was constructed and is shown in Figure 4c. The constructed total delays then were applied into DInSAR analyses. Figure 5 shows the deformation field derived from DInSAR analyses conducted with the error correction method described above.
temporal gap between MODIS water vapor and WRF simulation, the interpolated MODIS water vapor from the cloud free raw products and water vapor from WRF show similarities to some extents as shown in Figures 3b and 4b. Together with ZDD, based on the surface pressure in WRF and Equation (6), the total delay map was constructed and is shown in Figure 4c. The constructed total delays then were applied into DInSAR analyses. Figure 5 shows the deformation field derived from DInSAR analyses conducted with the error correction method described above.

Time Series Analysis with ENVISAT ASAR
StaMPS/MTI processor conducted the PS analysis with the master SAR image of 24 July 2009 and then SBAS analysis over three years period with the 18 images illustrated in Figure 6. The resulting velocity map and the mean precision dispersion maps derived from StaMPS/MTI processing are presented in Figures 7 and 8.

Time Series Analysis with ENVISAT ASAR
StaMPS/MTI processor conducted the PS analysis with the master SAR image of 24 July 2009 and then SBAS analysis over three years period with the 18 images illustrated in Figure 6. The resulting velocity map and the mean precision dispersion maps derived from StaMPS/MTI processing are presented in Figures 7 and 8.

Time Series Analysis with ENVISAT ASAR
StaMPS/MTI processor conducted the PS analysis with the master SAR image of 24 July 2009 and then SBAS analysis over three years period with the 18 images illustrated in Figure 6. The resulting velocity map and the mean precision dispersion maps derived from StaMPS/MTI processing are presented in Figures 7 and 8.  Table 2 for the image acquisition date.
For further interpretation, the area of timberline was extracted using the SAR classification scheme proposed by [62]. The boundaries are drawn using a red line in Figures 7b and 8b. It is observed that the mean velocity and its dispersion of scatterers over the area outside of timberline  are much higher than the stable scatterers within the timberlines. Therefore the following discussion was conducted based on the surface deformation detected inside the timberline area.   For further interpretation, the area of timberline was extracted using the SAR classification scheme proposed by [62]. The boundaries are drawn using a red line in Figures 7b and 8b. It is observed that the mean velocity and its dispersion of scatterers over the area outside of timberline are much higher than the stable scatterers within the timberlines. Therefore the following discussion was conducted based on the surface deformation detected inside the timberline area.

Discussion
Based on the DInSAR analysis with atmospheric correction and time-series processing, the errors and relevant issues are discussed below to achieve a reliable interpretation of volcanic activity.

DInSAR Performance
As shown in Figure 5, after the atmospheric correction, two pass DInSAR using ALOS PALSAR images demonstrated high deformation values. Particularly in mountain flank areas, the displacement rate was up to 15-20 cm/year, which is a very high value compared with the StaMPS/MTI output and close to the deformation observed at other activate volcanoes [63,64]. It was suspected that most of this deformation might be caused by the base DTM error remaining from the two pass DInSAR processing. To examine this error source, the data collected by Geoscience Laser Altimeter System (GLAS) on board the Ice, Cloud and land Elevation Satellite (ICESAT) was employed as the reference for inter-comparison. Due to the high vertical and

Discussion
Based on the DInSAR analysis with atmospheric correction and time-series processing, the errors and relevant issues are discussed below to achieve a reliable interpretation of volcanic activity.

DInSAR Performance
As shown in Figure 5, after the atmospheric correction, two pass DInSAR using ALOS PALSAR images demonstrated high deformation values. Particularly in mountain flank areas, the displacement rate was up to 15-20 cm/year, which is a very high value compared with the StaMPS/MTI output and close to the deformation observed at other activate volcanoes [63,64]. It was suspected that most of this deformation might be caused by the base DTM error remaining from the two pass DInSAR processing. To examine this error source, the data collected by Geoscience Laser Altimeter System (GLAS) on board the Ice, Cloud and land Elevation Satellite (ICESAT) was employed as the reference for inter-comparison. Due to the high vertical and planimetric accuracy (2 cm radially and 5 cm vertically [65], it is a highly reliable topographic product to assess the accuracy of the 3 arc-second resolution SRTM DEM used in the DInSAR processing over Mt. Baekdu. The comparison between SRTM DEM and ICESAT scanned profiles is shown in Figure 9. Along the two ICESAT track profiles, it was clearly observed that the horizontal and vertical geodetic offsets of STRM DTM existed with about 20 m height difference over the DInSAR processing area. Then the deformation error resulting from the residual height of the base DTM could be computed based on the equation: where Φ t is the phase difference via inaccurate base topography, λ is the wavelength of SAR sensor, B p is the perpendicular baseline, θ is the look angle, r is the range between sensor and target and z error is the base DTM error. As a result, the calculated deformation errors were up to 1.9 cm and 2.3 cm respectively in pairs 1 to 2 (about 10 cm/year, refer to Section 4.1 and Figure 5). Hence large portion of the deformation derived from the DInSAR processing originated from the geodetic offset of the base DTM.
where Фt is the phase difference via inaccurate base topography, λ is the wavelength of SAR sensor, Bp is the perpendicular baseline, θ is the look angle, r is the range between sensor and target and zerror is the base DTM error. As a result, the calculated deformation errors were up to 1.9 cm and 2.3 cm respectively in pairs 1 to 2 (about 10 cm/year, refer to Section 4.1 and Figure 5). Hence large portion of the deformation derived from the DInSAR processing originated from the geodetic offset of the base DTM. Together with the offset caused by the base DTM error, atmospheric error might be another unresolved error component. First, the discontinuities in error maps and Figure 5 should be noted. The discontinuities shown in the error maps ( Figure 5) in pair 1 and especially pair 2 clearly revealed the discontinuities around the summit of Mt. Baekdu, which might be due to the estimated frontal surfaces of WRF outputs. DInSAR deformation showed the unpredicted discontinuities by WRF in both pairs 1 and 2. These were artifacts derived from inaccurate WRF outputs as WRF has limited capability to forecast local weather instabilities. In addition, considering the locations over the boundaries of estimated discontinuities revealed in the WRF error map, deformation over some local area, for instance in the north eastern flank of pair 1 and the western flank in pair 2 (shown in Figure 5), were inferred as the result of a mismatch between WRF estimation and real atmosphere instability.
With all of the above unresolved phase delay errors, the two pass DInSAR results were unable to show precise genuine deformation patterns. There is only weak evidence of deformation on the eastern flank. Next, we use the StaMPS/MTI method to precisely track the deformation.

Interpretation of StaMPS/MTI Processing Results
In the StaMPS/MTI processing results, the obvious inflation on eastern flank was distinguished. Specifically over the site "a" indicated in Figure 7b, inflation up to 30 mm/year is shown. This deformation rate is significant and represent genuine surface change based on: (1) the coincidence with the deformation in the other SBAS campaign shown in Kim et al. [22]; (2) its Together with the offset caused by the base DTM error, atmospheric error might be another unresolved error component. First, the discontinuities in error maps and Figure 5 should be noted. The discontinuities shown in the error maps ( Figure 5) in pair 1 and especially pair 2 clearly revealed the discontinuities around the summit of Mt. Baekdu, which might be due to the estimated frontal surfaces of WRF outputs. DInSAR deformation showed the unpredicted discontinuities by WRF in both pairs 1 and 2. These were artifacts derived from inaccurate WRF outputs as WRF has limited capability to forecast local weather instabilities. In addition, considering the locations over the boundaries of estimated discontinuities revealed in the WRF error map, deformation over some local area, for instance in the north eastern flank of pair 1 and the western flank in pair 2 (shown in Figure 5), were inferred as the result of a mismatch between WRF estimation and real atmosphere instability.
With all of the above unresolved phase delay errors, the two pass DInSAR results were unable to show precise genuine deformation patterns. There is only weak evidence of deformation on the eastern flank. Next, we use the StaMPS/MTI method to precisely track the deformation.

Interpretation of StaMPS/MTI Processing Results
In the StaMPS/MTI processing results, the obvious inflation on eastern flank was distinguished. Specifically over the site "a" indicated in Figure 7b, inflation up to 30 mm/year is shown. This deformation rate is significant and represent genuine surface change based on: (1) the coincidence with the deformation in the other SBAS campaign shown in Kim et al. [22]; (2) its location in the non-vegetated area (hence this deformation was not the result of time varying scatterers); and (3) a relatively small velocity dispersion as shown in Figure 8. The sources of the other deformation in the small areas around the summit are not clearly identified as those are very localized pattern. In addition, the minor inflation on the site "b" is relatively obvious and reliable according to the low mean velocity dispersion shown in Figure 8.
The deformation over sites "c" on the southern flank and "d" o the western flank is not negligible; however it is not very clear in both the mean velocity and the velocity dispersion. The deformation of site "c" is more likely the noise of DInSAR processing or the deformation caused by other sources considering its localized distribution. Regarding the site "d", as the topographic slope in this area is higher than in any other place, the influence of decorrelation in steep topography may cause the erroneous measurement in phase difference as stated in [23]. The relatively low reliability in velocity dispersion was noticed over the western flank. However, since it was difficult to address the de-correlation problem originating from the vegetation canopy, such as shrub at high altitude, snow coverage or steep slopes which are dominant on the western flank, the high dispersion in this area was inferred as the effects of such natural environmental factors. Thus, together with the spatially small extents, the observed deformation in the western flank was not readily considered as real deformation signal. The site "a" is clearly distinguished in the well-defined uplift pattern as observed consistently in Figure 10. Although there are some fluctuations, the deformation over the site "a" is very clear for the entire observation period, which is also shown in Figure 11a. Moreover, the deformation on the site "b" on the northern flank area during the 2007-2008 period showed the small velocity dispersion in Figure 11b and is identified visually in Figure 10.
In contrast, it is found that the deformation detected over sites "c" on the southern flank and "d" on the western flank areas is ambiguous. As described by Cao et al. [66], the deformation extracted for site "c" might be local creep and it agrees with Figure 11c which shows slow subsidence with high velocity dispersion. As for site "d", the less obvious deformation with large deformation velocity dispersion shown in Figures 10 and 11d imply that the observed deformation over site "d" may not be actual surface uplift and subsidence but error in the DInSAR analysis.
Unfortunately, the reliable StaMPS/MTI observation areas in this study are not overlapped with [4] survey areas spatially and temporally especial ly "a" to "c" are located within North Korean territory; hence, the direct inter-comparison is not feasible. However, we can compare [4] survey work and the InSAR analyses conducted in this paper. Firstly, it was identified that there were still remained surface deformation in the 2007-2011 period after the long inflation during the 2002-2005 period observed in [4] ground survey, which was covered by the time series analysis using ENVISAT ASAR data and in the DInSAR analysis using ALOS PALSAR images. Specifically, the ENVISAT ASAR StaMPS/MTI analyses showed a clear evidence of surface deformation on the eastern flank of Mt. Baekdu, which was not covered by the on-site GPS and levelling survey conducted by [4] as the area is located within N. Korean territory. The phenomenon could also be observed in the SBAS analysis reported in [22]. All StaMPS/MTI observations identified quite obvious deformation (up to 15-30 mm) on the eastern flank. Thus, it ruled out the possibility of outliers in the DInSAR analyses originating from erroneous factors, especially considering the timberline extent and the small velocity dispersion. It should be noted that the circumscribed patterns of deformation connecting sites "a", "b" and "d" around the eastern summit, fit with the epicenters of swarms around the caldera lake [67]. These were described in Liu et al. (2001) [2]. Such circumscribed deformation patterns are also revealed in the atmospheric corrected two pass DInSAR results of pairs 1-3 shown in Figure 5. This implies that the deformation observed around the caldera lake is the consequences of ongoing geological unrest, not error components. It should be noted the levelling and GPS results conducted in [4] is not directly comparable to our StaMPS/MTI observations as the spatial extent, the temporal coverage and the measured direction were not within our test scope.
In order to infer the deformation source model of Mt. Baekdu and to validate the effectiveness of our DInSAR observations, we conducted model inversion incorporating the DInSAR observations for a specific time period. According to [5], the Earth Gravitational Model (EGM) 2008 [68] was applied which approximates a magma chamber along with several faults. Hence in this study we combined the Mogi [69] and Okada [70] models to simulate a comprehensive underground source causing surface deformation. The spherical Mogi model [69] handles circular surface displacement well, specifically when the source depth is much larger than the radius of the deformation cavity. Because of its simplicity it has been widely used for volcanic region inversion [71,72]. The rectangular Okada model [70] mimics fault behavior and is widely modeled in subduction and fault areas [73,74]. study we combined the Mogi [69] and Okada [70] models to simulate a comprehensive underground source causing surface deformation. The spherical Mogi model [69] handles circular surface displacement well, specifically when the source depth is much larger than the radius of the deformation cavity. Because of its simplicity it has been widely used for volcanic region inversion [71,72]. The rectangular Okada model [70] mimics fault behavior and is widely modeled in subduction and fault areas [73,74].   Both models were optimized using the Levenberg-Marquardt algorithm (LMA) which was implemented using Gauss-Newton iteration and loops in order to obtain a convergence of the gradient descent [75]. Initial model values for geometrical parameters are set to ranges published in previous studies [5]. For the Mogi model they include depth, projected planar location, and volume; for the Okada model they include, length, width, depth, degree of dip and rake, distance of strike and slip. Through the iterative calculation during the inversion, these parameters were optimized to minimize the differences between the model estimates and our DInSAR observations [76]. Once the processing was completed, parameters of the pressure source can be determined and described. Successful studies using inversion modelling have been published by, e.g., Yassen et al. [77] and Cheloni et al. [78].
We tested the inversions of all DInSAR observations. Since only a few StaMPS/MTI results showed clear surface deformation and produced the source parameters in reasonable ranges, we In the first scenario, a Mogi-type magma chamber of spherical shape in an elastic half-space [70] and four faults lines inferred from [5] were introduced. In this scenario the central magma chamber contributes to surface unrest by volume or pressure changes interacting with the faults lines as shown by, e.g., Feuillet et al. [79]. Figure 12 and Table 3 present the observed and modelled Both models were optimized using the Levenberg-Marquardt algorithm (LMA) which was implemented using Gauss-Newton iteration and loops in order to obtain a convergence of the gradient descent [75]. Initial model values for geometrical parameters are set to ranges published in previous studies [5]. For the Mogi model they include depth, projected planar location, and volume; for the Okada model they include, length, width, depth, degree of dip and rake, distance of strike and slip. Through the iterative calculation during the inversion, these parameters were optimized to minimize the differences between the model estimates and our DInSAR observations [76]. Once the processing was completed, parameters of the pressure source can be determined and described. Successful studies using inversion modelling have been published by, e.g., Yassen et al. [77] and Cheloni et al. [78].
We tested the inversions of all DInSAR observations. Since only a few StaMPS/MTI results showed clear surface deformation and produced the source parameters in reasonable ranges, we calculated inversions based on observed that occurred between 2007/08/  In the first scenario, a Mogi-type magma chamber of spherical shape in an elastic half-space [70] and four faults lines inferred from [5] were introduced. In this scenario the central magma chamber contributes to surface unrest by volume or pressure changes interacting with the faults lines as shown by, e.g., Feuillet et al. [79]. Figure 12 and Table 3 present the observed and modelled deformation and error residuals for this hypothesis. Although the modelled deformation values show significant differences when compared to the DInSAR observations, the spatial patterns between modelled and observed surface deformation are well comparable. To estimate the goodness of fit of the model more objectively, we introduced the statistical metric from Tizzani et al. [80] together with the conventional RMS error estimate between the observation and modelled displacement. Their statistical test is based on the experimental variogram model of deformation distribution. First, the normalized root mean square error (NRMSE) was calculated which provides a measure for the discrepancy between the observed and model variogram. It is defined as where obs and model are observational and model data and N is the number of data points. We also conducted a Chi-squared test between the observed and the model variogram and calculated χ 2 (see Tables 3 and 4). All cases show a comparable fitness although in the case of the 2007 observation it is slightly higher due to the possibility that the observed deformation possessed strong deviation. Volume change in the model magma chamber using Mogi source model parameters is also extracted using where ∆V is the volume change in the Mogi source, a is the source radius, ∆P is the pressure change in the source, z 0 is the depth of source and µ is the shear modulus [81] (see Table 3) .
Overall it appears that the deformation in the cases of the 2007 and 2008 observations is due to the inflation of the magma chamber and its corresponding interaction with fault lines. In contrast, the 2010 observations seem to show the depressurization of the magma chamber. Consequently the focal mechanism of faults lines ( Figure 12) display different patterns according to the magnitude and direction of volume changes of the Mogi source, in particular between the major inflation of the 2007 observation period and the deflation of the 2010 observation period. It is obvious that significant pressure change or any physical state migration of the magma chamber caused by, e.g., magma intrusion which can induce measureable and frequent major seismic activity, have ceased after 2005 [2]. Thus, the magma chamber did not significantly contribute to the observed surface deformation during the period of this study. Assuming the absence of major activity in the magma chamber, it was proposed that minor coulomb stress changes of the deep magma chamber induced the clamping and unclamping of fault planes as shown in the case studies by [79,82]. It is characterized by predominantly normal faulting in the focal mechanisms of the 2007 observation towards a tendency of strike-slip faulting in the 2008 observation. Deformation during the 2010 observation period is characterized by dominant strike-slip faulting which then might be directly connected to the deflation of the magma chamber ( Figure 12C).
It should be noted that the residuals between observed and modelled deformation shown in Figure 12 are significant although some deformation patterns such as discontinuities probably induced by faults can be roughly traceable in the observations. A possible explanation is the effect of phase errors induced by local atmospheric components such as local precipitation. However, in order to explain the distribution and migration of the residuals between observations and models and their temporal migrations, it is useful to infer the existence of additional small scale deformation sources. Such additional local faults might be trapdoor faults resulting from intense magmatic activity during the 2003-2005 period [2] as shown in Sierra Negra Volcano [83]. However, we found that an improvement of the model inversion including a magma chamber, major fault lines based on [5] and a number of unknown local deformation sources were not possible to obtain due to many unconstrained parameters uncertainties.  Instead, as the second scenario, we examined the model inversion including only locally distributed Okada type deformation sources assuming negligible contribution from a magma chamber and major fault lines as inferred from [5]. The assumption can be justified considering the decrease in seismic activity after 2005 [2] which implies a stabilization of the magma chamber. Therefore, the model inversion was conducted with only two rectangular dislocations in an elastic half space [70] and the initial parameter values inferred from the location of main area of deformation. The models incorporating the NW-SE and NE-SW trending fault lines provided the best results. Parameters derived from these models are shown in Table 4, together with observed and modelled deformation and error residual distributions as shown in Figure 13. Location of deformation sources in these models showed a better agreement with the depth of earthquake hypocenters (2-5 km) described in [4]. One hypothesis involving these local deformation models is that hydrothermal activity might have followed main inflation events and intensive seismic activity during the 2003-2005 period [2]. In such a scenario, close proximity between an aquifer and the magma chamber could have caused significant gas emission, including water vapor and induced relatively localized surface deformation along existing faults as observed in [84,85] or with the changes of pore pressure and rock temperature [86]. The focal mechanisms and the spatial patterns of deformation sources have been changing significantly across the three inversions and may also support a scenario of temporal migration of hydrothermally activated regions. The focal mechanisms clearly demonstrate strike-slip faulting and can be well discriminated from the cases of Figure 12A with predominant normal faulting and Figure 12B with a tendency towards strike-slip faulting. It seems that the error residual patterns in Figure 13 and higher error metrics in Table 4 (cf. to values in Table 3) require the introduction of more sophisticated deformation source models and mechanisms. However, the thermal anomaly around the summit observed in MODIS land temperature products [87] shown on Figure 14 strongly support the possibility of hydrothermally Instead, as the second scenario, we examined the model inversion including only locally distributed Okada type deformation sources assuming negligible contribution from a magma chamber and major fault lines as inferred from [5]. The assumption can be justified considering the decrease in seismic activity after 2005 [2] which implies a stabilization of the magma chamber. Therefore, the model inversion was conducted with only two rectangular dislocations in an elastic half space [70] and the initial parameter values inferred from the location of main area of deformation. The models incorporating the NW-SE and NE-SW trending fault lines provided the best results. Parameters derived from these models are shown in Table 4, together with observed and modelled deformation and error residual distributions as shown in Figure 13. Location of deformation sources in these models showed a better agreement with the depth of earthquake hypocenters (2-5 km) described in [4]. One hypothesis involving these local deformation models is that hydrothermal activity might have followed main inflation events and intensive seismic activity during the 2003-2005 period [2]. In such a scenario, close proximity between an aquifer and the magma chamber could have caused significant gas emission, including water vapor and induced relatively localized surface deformation along existing faults as observed in [84,85] or with the changes of pore pressure and rock temperature [86]. The focal mechanisms and the spatial patterns of deformation sources have been changing significantly across the three inversions and may also support a scenario of temporal migration of hydrothermally activated regions. The focal mechanisms clearly demonstrate strike-slip faulting and can be well discriminated from the cases of Figure 12A with predominant normal faulting and Figure 12B with a tendency towards strike-slip faulting. It seems that the error residual patterns in Figure 13 and higher error metrics in Table 4 (cf. to values in Table 3) require the introduction of more sophisticated deformation source models and mechanisms. However, the thermal anomaly around the summit observed in MODIS land temperature products [87] shown on Figure 14 strongly support the possibility of hydrothermally induced deformation. These data demonstrate a consistent thermal anomaly with up to 10 Celsius degree, in particular near the location of the south-eastern part of the caldera lake. We interpret these anomalies as possible outlets of hydrothermal activity based on the proximity to the modelled deformation sources.
Remote Sens. 2017, 9,138 19 of 25 induced deformation. These data demonstrate a consistent thermal anomaly with up to 10 Celsius degree, in particular near the location of the south-eastern part of the caldera lake. We interpret these anomalies as possible outlets of hydrothermal activity based on the proximity to the modelled deformation sources.

Conclusions and Future Work
Based on recent ground observations and seismic data, a number of indicators of potential volcanic activity around Mt. Baekdu could be identified, making a thorough re-evaluation of its potential volcanic risk critical. However, since a major part of Mt. Baekdu is not accessible for ground surveying and in situ studies, we focused on the implementation and validation of InSAR techniques that allowed us to monitor the target area both temporally and spatially.
To overcome the environmental conditions limiting the InSAR performance (such as potentially inaccurate base DTM, instability of atmospheric condition resulting in partial phase delay, dense vegetation canopy and steep slope), error elimination techniques, including two pass DInSAR with atmospheric correction and StaMPS/MTI methods were employed successively. As a result, it was verified that the local deformation identified in both approaches is not negligible. In induced deformation. These data demonstrate a consistent thermal anomaly with up to 10 Celsius degree, in particular near the location of the south-eastern part of the caldera lake. We interpret these anomalies as possible outlets of hydrothermal activity based on the proximity to the modelled deformation sources.

Conclusions and Future Work
Based on recent ground observations and seismic data, a number of indicators of potential volcanic activity around Mt. Baekdu could be identified, making a thorough re-evaluation of its potential volcanic risk critical. However, since a major part of Mt. Baekdu is not accessible for ground surveying and in situ studies, we focused on the implementation and validation of InSAR techniques that allowed us to monitor the target area both temporally and spatially.
To overcome the environmental conditions limiting the InSAR performance (such as potentially inaccurate base DTM, instability of atmospheric condition resulting in partial phase delay, dense vegetation canopy and steep slope), error elimination techniques, including two pass DInSAR with atmospheric correction and StaMPS/MTI methods were employed successively. As a result, it was verified that the local deformation identified in both approaches is not negligible. In

Conclusions and Future Work
Based on recent ground observations and seismic data, a number of indicators of potential volcanic activity around Mt. Baekdu could be identified, making a thorough re-evaluation of its potential volcanic risk critical. However, since a major part of Mt. Baekdu is not accessible for ground surveying and in situ studies, we focused on the implementation and validation of InSAR techniques that allowed us to monitor the target area both temporally and spatially.
To overcome the environmental conditions limiting the InSAR performance (such as potentially inaccurate base DTM, instability of atmospheric condition resulting in partial phase delay, dense vegetation canopy and steep slope), error elimination techniques, including two pass DInSAR with atmospheric correction and StaMPS/MTI methods were employed successively. As a result, it was verified that the local deformation identified in both approaches is not negligible. In particular, the hybrid interferogram stacking algorithms developed by [34] revealed significant surface inflation up to 20 mm/year supposedly around faults over the caldera lakes. To explain this deformation, we propose two hypotheses: (1) faults interacting with the magma chamber and (2) hydrothermal activity and effects around local faults in the flank over Mt. Baekdu. Although both hypotheses do not fully explain the complicated deformation pattern identified by DInSAR observations, the temporal migration of the modelled deformation source, calculated focal mechanism and thermal anomaly detected by MODIS land products support a scenario of recently ongoing volcanic activity. A unified scenario combining both mechanisms, i.e., vestigated but the final model inversion is not yet determined due to the difficulties in building more complex models, in particular with remaining atmospheric errors and insufficient DInSAR observations. Once reliable long-term observations are derived from ground observation and also remote sensing analyses, an identification and a detailed model of the deformation source underneath Mt. Baekdu could be derived. To achieve this purpose, the approaches combining InSAR analyses from different observation conditions and techniques [88,89] and/or data fusion with GPS monitoring [90], are suggested to be an appropriate bases for complex model inversion to infer various deformation sources.
In this work the effectiveness of DInSAR time series analysis over Mt. Baekdu has been demonstrated and exploitation of all possible space-borne SAR missions to determine the model that best fit the data over Mt. Baekdu is necessary. However, advanced error correction methods need be developed and introduced to improve the results. Although time series using StaMPS/MTI with C-band ENVISAT ASAR data produced relatively reliable deformation measurements, it is suggested to introduce long wavelength L-band DInSAR time series for highly intensive measurements over vegetated and/or snow covered topography at Mt. Baekdu. In addition, the base topography error issue should also be considered when time series analysis with L-band ALOS PALSAR is employed due to its temporal dependency on the perpendicular baseline [91]. In order to address these issues, approaches proposed by [91] and [92] can be introduced to make base DTM corrections for constant monitoring over Mt. Baekdu. Although it was shown that two-pass DInSAR compensated with external error estimation by weather forecasting models is not very effective in such severe atmospheric turbulence cases, it is still worthwhile to employ this method to investigate relatively short term deformation. The higher resolution and the improved accuracy of the WRF model to forecast the local climate instability over mountain areas such as the approach described in [93] is essential in that case.
The replacement of space-borne SAR sensors from ERS, ENVISAT ASAR and ALOS PALSAR to C-band Sentinel-1 of ESA, X-band high resolution SAR, and L-band ALOS PALSAR-2 provide enhanced opportunities to monitor Mt. Baekdu with improved spatial and temporal resolutions at various wavelengths. Together with the new SAR missions and improved DInSAR techniques described in this work, further monitoring over Mt. Baekdu could be continued and results described in this study, in particular with respect to various atmospheric error components, will be of great benefit.