Time Series Surface Deformation of Changbaishan Volcano Based on Sentinel-1B SAR Data and Its Geological Significance

Abstract

Monitoring the surface deformation is of great significance, in order to explore the activity and geophysical features of the underground deep pressure source in the volcanic regions. In this study, the time series surface deformation of the Changbaishan volcano is retrieved via Sentinel-1B SAR data, using the SBAS-InSAR method. The main results are as follows. (1) The mean surface deformation velocity in the Changbaishan volcano is uplifted as a whole, while the uplift is locally distributed, which shows a strong correlation with faults. (2) The time series surface deformation of the Changbaishan volcano indicates an apparently seasonal change. (3) The cumulative surface deformation shows a strong correlation with the maximal magnitude and number of annual earthquakes, and it is likely dominated by the maximal magnitude of the annual earthquakes. (4) The single Mogi source model is appropriate to evaluate the deep pressure source in the Changbaishan volcano, constrained by the calculated surface deformation. The optimal estimated depth of the magma chamber is about 6.2 km, and the volume is increased by about 3.2 × 10⁶ m³. According to the time series surface deformation, it is concluded that the tectonic activity and faults, related to the deep pressure source, are pretty active in the Changbaishan volcano.

1. Introduction

The Changbaishan Tianchi volcano (Changbaishan volcano) is located on the border between China and North Korea, the only large-scale Cenozoic multi-genetic central composite stratovolcano. It has experienced many violent eruptions in history, and the eruption in 946 AD is considered one of the largest volcanic eruptions in the past two thousand years, which is still thought as the most potentially hazardous active volcano in China [1,2,3,4,5].

Currently, the magmatic activity and physical properties of the magma chamber beneath the Changbaishan volcano are mainly studied with geophysical methods, including seismic, gravity, and geomagnetic exploration. The seismic tomography is used to study the origin of deep mantle volcanoes, the seismic prospecting method is employed to invert the crustal velocity structure, and the magneto-telluric sounding method is undertaken to invert the crustal resistivity structure; all these methods prove the existence of an active magma chamber beneath the Changbaishan volcano, formed by the subduction and dehydration of the Pacific plate [6,7,8,9,10,11,12,13,14,15,16]. The Bouguer gravity anomalies imply that there occurs a local mass loss, with low values centered at the Changbaishan volcano, indicating the existence of a magma chamber in the deep crust [8]. The in-situ measurements from 2002 to 2005 show that there is a strong correlation between the seismic intensity and magmatic activity in the Changbaishan area [17,18,19,20,21,22]. Thus, it is necessary to enhance the observation of the magmatic activity in the Changbaishan volcano.

In addition to geophysical exploration methods, remote sensing technology is also an important technique of volcanic activity monitoring. A widely used remote sensing sensor, synthetic aperture radar (SAR), is widely used not only in SAR change detection [23,24], but also in SAR modeling and SAR deformation monitoring of volcanoes [25,26]. The magmatic activity will bring about surface deformation in the volcanic regions [27,28,29]. Thus, besides the geophysical methods, the surface deformation is one of the important ways to indirectly monitor the magmatic activity of the magma chamber [30,31,32,33,34,35]. Since 2002, GPS and leveling surveys have been used to monitor the surface deformation in the Changbaishan volcano [36]. However, the operations of the field observation are critically limited by the complex surface topography and severe climate in the Changbaishan volcano. Comparatively, an interferometric synthetic aperture radar (InSAR) technique has commonly been considered one of the most profitable and effective ways to monitor the surface deformation in volcanic regions [37,38,39,40], which have also been used to study the Changbaishan volcano [30,31,32,33,34,35] (

Table 1. Current status in studying surface deformation in the Changbaishan volcano using InSAR method.

	Period	Dataset	Quantity of Datasets	Method	Maximum Surface Deformation	Pressure Source Inversion
Kim and Won [41]	1992–1998	ERS-1/2 and JERS-1	10/41	DInSAR	3 mm/year	-
Chen et al. [30]	1995–1998	JERS-1	2	DInSAR	almost no deformation	Double Mogi source
Han et al. [31]	1994–1998	JERS-1	4	DInSAR	5 mm/year	-
Tang P. et al. [35]	1992–1998	JERS-1	19	PS-InSAR	6 mm/year	-
	2007–2010	ENVISAT ASAR	14	PS-InSAR	3 mm/year	-
Kim et al. [34]	2007–2010	ENVISAT ASAR/ALOS PALSAR	19/4	StaMPS/DInSAR	30 mm/year	Single Mogi source/Okada dislocation
He et al. [32]	2006–2011	ALOS PALSAR	24	PS-InSAR	10 mm/year	Single Mogi source
Trasatti et al. [42]	2018–2020	ALOS-2	19	PS-InSAR	20 mm/year	Prolate spheroid

).

With the JERS-1 and ERS imagery, Kim and Won [41] firstly evaluated the surface information features in the Changbaishan volcano from 1992 to 1998, which indicates an average 3 mm/year uplift. Han et al. [31] assessed that the annual surface deformation of the Changbaishan volcano is about 5 mm/year, from 1994 to 1998, with JERS images, while Tang P. et al. [35] proposed a relatively larger annual surface deformation in this period, about 6 mm/year, using the permanent scattering InSAR (PS-InSAR) method. Unlike others, Chen et al. [30] found that the Changbaishan volcano had almost no deformation from 1995 to 1998. Tang P. et al. [35] obtained an average 3 mm/year uplift, from 2007 to 2010. Mainly in this period, He et al. [32] found that the uplift is about 10 mm/year using 24 ALOS PALSAR images, and Kim et al. [34] proposed a clear uplift up to 30 mm/year in the eastern slope of the Changbaishan volcano with ALOS PALSAR data and ENVISAT ASAR data. Using 19 ALOS-2 images from ascending orbit, provided by the Japanese Aerospatial Exploration Agency, Trasatti et al. [42] detected a 20 mm/year movement elongated in NW-SE direction on the southeastern flank and a −20 mm/year movement on the southwestern flank of the Changbaishan volcano from 2018 to 2020. Additionally, the depth and geometric properties of the deep pressure source are evaluated with the reckoned surface deformation by He et al. [32] and Kim et al. [34]. These studies bring a fresh view of the volcanic activities in the Changbaishan volcano.

However, the evaluated surface deformations are apparently different by several researchers, which is mainly brought by the following causes. First, the atmospheric error is hard to reduce in the traditional methods for the great elevation change and complex climate in the Changbaishan volcano. Second, the orbital periods of the used JERS-1, ASAR, and ALOS PALSAR sensors are different, and a stable data source with an excellent revisiting period is indispensable to assess the surface deformation features in the Changbaishan volcano.

As a new data source, Sentinel-1B SAR data have the good spatial resolution, extremely high orbit determination accuracy, and shorter revisit period, which will inevitably improve the monitoring of the surface deformation in the Changbaishan volcano. Additionally, the small baseline subset InSAR (SBAS-InSAR) method is feasible to significantly reduce the atmospheric error by the free combination of limited short time baseline and spatial baseline [43]. Therefore, the surface deformation in the Changbaishan volcano is evaluated with the SBAS-InSAR method, using the Sentinel-1B SAR data.

2. Study Area and Data

2.1. Regional Geology

The Changbaishan volcano is located on the border between China and North Korea in the east of Jilin Province, ranging from 41°30′N to 42°36′N and 127°22′E to 128°41′E (Figure 1). The highest altitude is about 2749.2 m, at Mount Baitou [2,3]. Here, Wangtian’e volcano is located approximately 30 km southwest, and Nampotae volcano is about 40 km northwest of the Changbaishan volcano. The three volcanoes belong to the evolutionary sequence of the Changbaishan volcano group, and the youngest active one is the Changbaishan volcano. It has a temperate continental mountain climate, with an average annual rainfall of 1333 mm. The snow cover always begins in October and ends in April (of the following year). More than 85% of the study area is covered by a dense forest, and there is a typical division of the vertical vegetation zones.

The Changbaishan volcano is located to the east of the Dunhua-Mishan fault in the outermost part of the Northeast Rift System, the front of the subduction zone of the Western Pacific Plate, and the northern edge of the ancient Sino- Korean craton fault [2,44,45]. There are mainly ten faults marked f1 to f10, as shown in Figure 1 [15,45,46,47,48]. Fault f1 is originated from Baishanzhen to Kimchaek in NW trending, which is nearly upright and 320 km in length. Fault f2 is from Liudaogou through the Changbaishan volcano until Zhenfengshan in NE trending, with a distance of 150 km, which is a compressive-shear fault. Faults f1 and f2 are the most significant geologic structures in the Changbaishan region, which dominated the shape of the Changbaishan volcano. Fault f3, a strike-slip fault with a length of over 120 km, is extended from Mayihe to Liangjiang in NE trending. Fault f4 is from Manjiang to Baishan in NWW trending, and a fault cliff in V shape is developed along this fault, forming the main valley of the Manjiang River. Fault f5, with a length of about 50 km, is from Songshan to Tianchi in NWW trending, a radial fault in the west of the Changbaishan volcano. Fault f6 is a tensile fault with NNE trending, which is from Tianchi to Hongtoushanzi. Faults f7 is from Tianchi to Tumen River, spreading along the NW direction, controlling the Pleistocene basalt volcanic group in the study area. Faults f8, f9, and f10 are interpreted with the remote sensing technique by Chen [46].

2.2. Sentinel-1B SAR Data

The Sentinel-1B satellite was launched on 25 April 2016, and it began to acquire SAR data products in C-band at a 5.4 GHz frequency in August 2016 by the European Space Agency (ESA) Copernicus program. As an earth observation satellite, the 12-day revisit period is reasonably short, and the orbit determination ephemeris is highly accurate, making Sentinel-1B SAR data widely used to monitor the surface deformation in extensive regions [49,50].

To obtain the time series surface deformation of the Changbaishan volcano, 75 descending orbit images were used continuously from 2016 to 2020. The detailed information of the selected data is listed in Appendix A 

Table A1. Available Sentinel-1B data of the study area, from 2016 to 2020.

No.	Acquisition Time	Absolute Orbit	Path	Perpendicular Baseline (m)	No.	Acquisition Time	Absolute Orbit	Path	Perpendicular Baseline (m)
0	26 November 2016	003135	134	93.8761	38	24 August 2018	012410	134	49.3319
1	8 December 2016	003310	134	74.4027	39	5 September 2018	012585	134	53.1985
2	20 December 2016	003485	134	93.3301	40	29 September 2018	012935	134	49.6745
3	1 January 2017	003360	134	144.857	41	11 October 2018	013110	134	134.333
4	25 January 2017	004010	134	106.654	42	23 October 2018	013285	134	144.707
5	18 February 2017	004360	134	42.798	43	4 November 2018	013460	134	91.89
6	2 March 2017	004535	134	−24.2414	44	16 November 2018	013635	134	78.573
7	14 March 2017	004710	134	107.537	45	10 December 2018	013985	134	101.761
8	26 March 2017	004885	134	110.998	46	22 December 2018	014160	134	174.279
9	7 April 2017	005060	134	66.7698	47	15 January 2019	014510	134	354.1653
10	19 April 2017	005235	134	35.2373	48	8 February 2019	014860	134	119.087
11	1 May 2017	005410	134	57.174	49	20 February 2019	01 5035	134	153.503
12	6 June 2017	005935	134	99.3488	50	4 March 2019	015210	134	138.687
13	18 June 2017	006110	134	24.6913	51	16 March 2019	015385	134	9.10079
14	30 June 2017	006285	134	−11.1132	52	9 April 2019	015735	134	44.1878
15	24 July 2017	006635	134	56.8439	53	21 April 2019	015910	134	98.6213
16	5 August 2017	006810	134	110.937	54	3 May 2019	016085	134	167.254
17	17 August 2017	006985	134	−11.6751	55	27 May 2019	016435	134	68.9393
18	29 August 2017	007160	134	24.4308	56	8 June 2019	016610	134	−14.0697
19	10 September 2017	007335	134	55.7565	57	20 June 2019	016785	134	60.5642
20	22 September 2017	007510	134	112.695	58	14 July 2019	017135	134	97.993
21	4 October 2017	007685	134	87.7212	59	26 July 2019	017310	134	25.8509
22	28 October 2017	008035	134	69.9597	60	7 August 2019	017485	134	12.6419
23	21 November 2017	008385	134	106.174	61	31 August 2019	017835	134	72.3225
24	15 December 2017	008735	134	133.51	62	12 September 2019	018010	134	77.0498
25	8 January 2018	009085	134	79.7278	63	24 September 2019	018185	134	86.2929
26	1 February 2018	009435	134	145.222	64	6 October 2019	018360	134	27.5584
27	13 February 2018	009610	134	147.873	65	18 October 2019	018535	134	77.5763
28	25 February 2018	009785	134	24.2284	66	11 November 2019	018885	134	114.73
29	21 March 2018	010135	134	−57.343	67	17 December 2019	019410	134	86.3005
30	2 April 2018	010310	134	104.749	68	29 December 2019	019585	134	142.613
31	8 May 2018	010835	134	17.6941	69	10 January 2020	019760	134	161.581
32	20 May 2018	011010	134	−10.57	70	22 January 2020	019935	134	125.07
33	1 June 2018	011185	134	−24.3215	71	3 February 2020	020110	134	79.863
34	25 June 2018	011535	134	69.2929	72	15 February 2020	020285	134	37.91
35	7 July 2018	011710	134	40.1216	73	10 March 2020	020635	134	138.783
36	31 July 2018	012060	134	0	74	22 March 2020	020810	134	121.937
37	12 August 2018	012235	134	53.4225	-	-	-	-	-

. All the data are single look complex (SLC) products, with descending orbit (path = 134), interferometric wide swath (IW) mode, and VV polarization. The spatial resolution of the data is 20 m × 5 m. Moreover, the time interval of continuously acquired images is limited to one revisit period (12 days) to three revisit periods (36 days) to maintain a high coherence. Then, with the available 75 scenes of the SAR data, the perpendicular baseline is less than 180 m, with an average of 83 m, according to Appendix A Table A1. The above parameters are better compared with the SAR data used by Tang P. et al. [35] and Kim et al. [33,34], when studying the surface deformation in the Changbaishan volcano.

Additionally, ALOS World 3D-30 m (AW3D30) data, global digital surface model (DSM) data released by Japan Aerospace Exploration Agency (JAXA) in May 2015, have been used to coregister the SAR images accurately and to correct interferograms brought by the topography. The spatial resolution of the data is 30 m, and the elevation accuracy is 5 m [51]. To avoid the baseline error, caused by the error-containing orbit information, the POD precise orbit ephemerides data are also used, whose orbit determination accuracy is better than 5 cm [51].

3. Data Processing

3.1. Subsection

SBAS-InSAR is systematically and widely used to retrieve the deformation field of the volcanic region, and it can document centimeter- to millimeter-scale slow surface deformation [52]. By the free combination of 75 SAR images, with short time intervals and perpendicular baselines, the influences of spatiotemporal decoherence and atmospheric delay can be well-suppressed with the SBAS-InSAR technique.

In this study, the SBAS-InSAR method was performed with SARSCAPE software, which was developed by the SARMAP AG Company. At first, by selecting a definite temporal-spatial baseline threshold (according to Figure 2, the time baseline is 430 days and perpendicular baseline is 180 m) and considering the Doppler center frequency difference, the data set is divided into different subsets by Delaunay triangulation to generate 1533 interferometric pairs (Figure 2a). Then, the accuracy of the sub-pixel level was achieved by using the external DEM data for the auxiliary accurate co-registration. Additionally, each pair of matched images was imaged by interference, and interferogram and coherence coefficient maps were obtained. Finally, after evaluating the quality of all interferograms and coherence coefficient maps, the interferometric pairs with low coherence were removed, and only 309 interferometric pairs were kept (Figure 2b).

The reference DEM was used to remove the flat effect phase and topographic phase, and an adaptive power spectrum filter [53] was adopted to reduce the phase noise in the interferograms before phase unwrapping with the Delaunay minimum cost flow (Delaunay MCF) [54] method.

A total of 86 ground control points were selected in the regions with high coherence coefficients, which were used to estimate and remove the phase changes caused by orbit error, the residual topography phase, and slope phase in unwrapping phase superposition, using the cubic polynomial refinement method. Then, a linear model between deformation velocity and elevation coefficient, based on highly coherent pixels, was established to derive the residual height and preliminary deformation, and the singular value decomposition method was applied to search the least-squares solution for each coherent pixel, as well as to estimate the nonlinear deformation. After eliminating the influence of atmospheric delay through a spatial-temporal filter, the cumulative surface deformation (CSD) and mean deformation velocity (MDV) along the line of sight (LOS) were obtained, at last (Figure 3).

3.2. Mean Deformation Velocity (MDV) in Changbaishan Volcano

Figure 3 shows that the number of the points with the estimated MDV are densely distributed on the northern and eastern slopes of the Changbaishan volcano, which is brought by the observation mode of the Sentinel-1B instrument orientating in this direction.

However, Figure 3 shows that the MDV is heavily dispersed in the Changbaishan volcano. Generally, the relatively lower absolute values mainly occur in the regions surrounding Tianchi caldera, while the relatively larger values (absolute values) occur in the region far from the caldera. Particularly, there also exist abnormal values in the caldera lake due to the specular reflection of SAR signals by water surface. Additionally, in the regions relatively far from Tianchi caldera, the abnormally large or small deformation values are irregularly distributed. Considering that the ground surface in these regions is densely covered by the vegetation and the climate in the Changbaishan volcano is rather complex, the influence of the atmospheric effect and surface topography will bring a severely temporal de-coherence of the scatterers and result in the dispersity of the surface deformation results. Thus, it is indispensable to refine the reliability of the calculated MDV.

3.3. Reliability Analysis

Here, the standard deviation of coherence is introduced to refine the MDV obtained with SBAS-InSAR method, which can be used to evaluate the quality of the retrieved results [33,55]. The corresponding formula is expressed as:

$$Pr{e}_v=\frac{\lambda }{4\pi }\sqrt{\frac{1-{\gamma }^2}{2{\gamma }^2}}$$

(1)

where λ is the SAR wavelength, γ is interferometric coherence, and Pre_v is the estimated value of the MDV accuracy. Based on Equation (1), the Pre_v of every inverted MDV point in the Changbaishan volcano can be calculated, as is shown in Figure 4a.

Figure 4 shows the relatively lower Pre_v values with green color and extending towards the southeast direction surrounding Tianchi caldera and the rather large values in the remaining place largely with red and yellow colors. That is, the MDV around Tianchi caldera are relatively higher accuracy.

To identify the surface deformation with high accuracy, the statistical histogram of the Pre_v is generated (Figure 4b), which shows a bimodal characteristic. The first mode ranges from 0.34 to about 1.176 mm, with the peak at 0.64 mm, representing the deformation values with higher accuracy. The second mode is from 1.176 to 3.08 mm with the 2.07 mm peak, representing the deformation values with poor accuracy. Thus, the surface deformation map of the Changbaishan volcano is generated with the first-mode Pre_v constraint (Figure 5).

Figure 5 shows that, after Pre_v constraint, the abnormal deformation values in the caldera lake and the regions relatively far from Tianchi caldera are totally eliminated in Figure 5. Compared to the previous studies [31,32,33,34,35], several new phenomena should be mentioned.

First, the uplift is different in dynamic ranges. Our result shows that the uplift velocity is generally from 0 to 5 mm/year. The maximum uplift velocity was 14.7 mm/year, which only existed in some isolated spots on the western slope. Comparatively, the uplift is generally from 0 to 35 mm by Kim et al. [33], 0 yearto 30 mm/year by Kim et al. [34], 0 yearto 5 mm/year by Han et al. [31], 0 yearto 6 mm/year by Tang P. et al. [35], 0 yearto 10 mm/year by He et al. [32], and 0 yearto 20 mm/year by Trasatti et al. [42]. Generally, the dynamic range of the uplift in this study is clearly smaller than that obtained by Kim et al. [33,34], but similar to that estimated by the others.

Second, Kim et al. [33,34] found a region with clearly uplift in the eastern slope, marked by A centered at (41°58′N, 128°08′E), which is clearly expressed in our study. Kim et al. [33] proposed that the uplift is about 30 mm/year. Here, Trasatti et al. [42] provided an up to 15 mm/year uplift, while it is only about 3 mm/year in our results. Besides Region A, more regions with similar uplift are found in this study, and they mainly exist in Regions B (42°1′53″N, 128°4′53″E), C (42°0′51″N, 128°2′16″E), and D (41°58′0″N, 128°4′10″E).

Third, besides the uplift, there also occur subsidence in abundant regions, which is seldom mentioned in the previous studies. The most apparent subsidence occurs in the eastern inner wall of Tianchi caldera (Region F, centered at (42°0′2″N, 128°4′39″E)), with obvious landslide characteristics. The area is about 4.45 km², and the maximal MDV is nearly −11.3 mm/year. Additionally, subsidence occurs in multiple regions with considerable area in the eastern slope of the Changbaishan volcano and a few patches with relatively small areas in the southern and northern slopes, including Region E (41°59′44″N, 128°6′38″E).

The reason to these phenomena will be further discussed in the following section.

4. Discussion

4.1. Spatial Distributions of Surface Deformation Fields

Generally, Figure 5 shows a fresh view of the surface deformation in the Changbaishan volcano, which is expressed in the following aspects.

(1)

The uplift is locally distributed.

Figure 5 shows that there occur abundant regions with uplift. Besides the uplift mentioned by Kim et al. [33,34] in Region A, the uplift mainly occurs in Regions B (northern slope) and C (western slope) of the Changbaishan volcano, and the deformation velocities of which are almost 2.2 mm/year, larger than that in Region A. Nearly all the most apparent uplifts occur in Region C in the western slope, and the mean uplift velocity is 4.37 mm/year. The highest uplift is up to 13.1 mm/year in Spots d (128°3′4″E, 42°1′46″N) and e (128°1′59″E, 41°59′29″N). Moreover, Region D in the southeastern slope is also uplifted, and its MDV is about 2.4 mm/year. In short, the uplifted extent is high in the east part of the study area in Region A, relatively low in the middle region located in the eastern and southern slopes, but highest in the western and northern slopes.

Interestingly, we notice that the boundaries of the locally distributed uplifts are related to the faults. Bounded by fault f2, the uplift degree of the northwestern part was generally higher than that of the southeastern part. Besides, along fault f1, there occur several patches with almost highest uplift in the northern part, while the east boundary of Region D and west boundary of Region A approximate the fault. Additionally, the surface deformation is clearly different between the east and west sides of fault f6.

To better understand the relationship between the surface deformation and the influence of the tectonic activities, the CSDs in multiple time spans are generated, as shown in Figure 6. Figure 6 verifies the influence of the tectonic activities on the evolution of the surface deformation, which is expressed as follows.

First, fault f2 is essential in controlling the variation of the locally distributed uplift. In December 2016 and January 2017, the difference is unclear between the two sides of the fault. Whereas from June 2017, the CSD in the west side is clearly higher than that in the east side, and the difference becomes fairly apparent in February 2020. Additionally, in February 2020, the uplift in the west side is almost the highest in the Changbaishan volcano, indicating that fault f2 plays a decisive role in the surface deformation in the study area. This phenomenon also agrees well that fault f2 is a significant dip-slip reverse fault in the Changbaishan volcano [45,46,47,48].

Second, fault f6 should also be in active. Bounded by fault f6, abundant spots with considerably large surface deformation occur in the west side and the local moderate uplift occurs in the east side, which is clearly expressed by the CSDs from June 2017. This means the strong correlation between fault f6 the locally distributed CSDs in the Changbaishan volcano, and its importance is just next to fault f2. However, the importance of fault f6 is not mentioned before, indicating that this should be an ignored significant fault.

Finally, the relationship between fault f1 and the locally distributed CSDs is not clear. Jin and Zhang [47], Liu et al. [48], Li et al. [45], Trasatti et al. [42], and Chen [46] proposed that fault f1 is one of the crucial faults in the Changbaishan volcano. Figure 6 shows that the difference of the CSDs along fault f1 is not clear. However, we mentioned that there occur abundant spots with relatively highest CSD in the northwestern section of fault f1, and the eastern section passes through the northeast boundary of Region D and southwestern boundary of Region A. Additionally, Using 19 ALOS-2 images from ascending orbit, Trasatti et al. [42] proposed that the fault is most active from 2018 to 2020. Here, considering that the used Sentinel-1B SAR data are only along a descending orbit, the deformation pattern retrieved could be influenced by the LOS in this study. Due to the parallel between the LOS vector and the fault strike, the locally distributed CSDs along fault f1 is not clear with respect to the other faults. Therefore, we ranked fault f1 as the third most active fault.

Therefore, the strong correlation between the surface deformation and the faults means that the tectonic activities play an essential role in controlling the current surface deformation of the Changbaishan volcano. Additionally, from the deformation quantity as expressed in Figure 5 and Figure 6, it can be concluded that fault f2 plays a controlling role on the surface deformation in the Changbaishan volcano, followed by faults f6 and f1.

(2)

There occurs abundant regions with subsidence.

Apparently, the subsidence is widely distributed in the study area. On a local scale, it is mainly distributed in the eastern slope of the study area, including Regions E and F, and the largest subsidence is up to −11.3 mm/year. In an isolated pattern, it is commonly distributed in a linear shape along the valleys. Particularly in Region E, the subsidence in local scale is extended in the linear shape along the river gully. Thus, erosion should be the main possibility leading to subsidence in the Changbaishan volcano.

Li et al. [45,56] and Liu et al. [2] suggested that the Changbaishan volcano is covered with pumice pyroclastic rocks of the Yuanchi and Tianwen periods, as well as the alkali flow pyroclastic rocks of the Qixiangzhan period, which were easy to be eroded and transported again under the condition of water saturation and caused volcanic debris flow and other disasters [57]. Thus, the strong correlation between the shape and distribution of the subsidence and lithology of the strata verifies that the subsidence is mainly brought by the erosion of the strata, which is not related to the activity of the deep pressure source in the Changbaishan volcano.

Additionally, on a local scale, the subsidence is widely distributed in the eastern slope. Comparatively, the uplift occurs in its east-most region and the western slope of the Changbaishan volcano. The surface deformation of such kind is probably related to the relative movements between Asia and Pacific Ocean plates, which deserves to be further studied in a more extensive scale.

4.2. Temporal Changes of Surface Deformation Fields

To further understand the time series characteristics of surface deformation in the study area, nine spots representing typical locations with evident uplift and subsidences in the Changbaishan volcano (Figure 5, Spots a–i) are selected to make the variation curves of CSD with time (Figure 7). Additionally, the first- (Y₁) and third-order (Y₂) polynomial fitting schemes of the CSDs with time are performed as an auxiliary to understand the time series characteristics of the surface deformation. The first-order polynomial fitting curve clearly shows the overall trend of the CSDs, and the third-order polynomial fitting curve presents detailed information of the CSDs. Several meaningful findings are presented, as follows.

(1)

The time series cumulative surface deformation is seasonal.

Figure 7 shows that the CSDs of Spots c, f, g, and h indicate a similar regularity. From November 2016 to August 2017, they all showed uplift, and the amount was less than 20 mm. Then, by March 2018, the surface deformation was the subsidence trend. The minimum subsidence was 2.27 mm at Spot g, while the maximal subsidence was only 7.36 mm at Spot h. From March 2018 to August 2018, they all showed the uplift, with a lifting amount of more than 6.09 mm. From September 2018 to March 2019, they indicate the subsidence trend, ranging from 3.12 mm to 8.71 mm. After that, Spots f, g, and h show an obvious uplift trend. Though Spot c showed an interval peak in August 2019, it presents an uplift trend on the whole.

Generally, the time series CSDs are apparently seasonal. The subsidence mainly occurs from October to February or March of the following year, while the uplift occurs from April to August or September.

This seasonal variation also occurs in the specific years in Spots d and e with obvious uplift. Additionally, though the time series CSDs in Spots a and b share apparent subsidence trends, the local CSD peaks primarily occur from April to August.

Regarding that the snow and ice cover of the mountain is seasonal, it is probable the existence of the snow and ice that brought the seasonal change of the time series CSDs. Whereas, Figure 7 indicates that, even if there exists the seasonal change of the time series CSDs, its quantity is slight compared to the total CSDs from 2016 to 2020. Thus, the retrieved time series CSDs are effective in this study.

(2)

The time series cumulative surface deformation is probably correlated with the tectonic activity.

Figure 7 shows that the variation of CSD with time probably reflects the tectonic activity of the Changbaishan volcano to some extent. From the overall performances (first-order fitting curves) and third-order polynomial fitting curves of the times series CSD, the Changbaishan volcano should be active from October 2016 to February 2017, and the CSD has increased significantly in this period. Thereafter, it should be quiet until May 2019, and the overall uplift of each point is not obvious. There is a slight downward trend at Spots d and g due to seasonal changes and stratum erosion. After May 2019, there probably occurs a new active period, and all spots show an apparent uplift trend.

To pursue the relationship between time series CSD and the tectonic activity, we superimposed the maximal magnitude and the number of annual earthquakes on Figure 7d,e,g,i with time series CSDs in similar pattern. Compared to the time series CSD and the maximal magnitude and the number of annual earthquakes, a meaningful finding is that the time series CSD is likely related to the earthquakes that happened in the Changbaishan volcano. Here, the annual maximal magnitude of the earthquakes increases from M_L 1.2 in 2016 to M_L 2.2 in 2018, and the number increases from 4 times in 2016 to 49 times in 2018; correspondingly, the maximal uplift of the fitted CSD has a 31 mm increase in this period. After that, the maximal magnitude of annual earthquakes decreases to M_L 1.9, and the number decreases to 19 times in 2019; simultaneously, the fitted CSD declined to largely 2 mm in this period. Then, from 2019 to 2020, the maximal magnitude keeps almost constant, but the number increases from 19 times to 99 times, while the maximal increase of the CSD is about 13 mm.

Thus, the time series CSD of the Changbaishan volcano is probably correlated with the maximal magnitude and the number of annual earthquakes. As mentioned above, from 2016 to 2019, the change of the fitted CSD is positively related to the maximal magnitude and the number of annual earthquakes, though the correlation is not apparent in the year from 2019 to 2020.

To further understand the relationship of the CSD and the earthquakes, the maximal magnitude and the number of annual earthquakes from 1999 to 2020, monitored by Changbaishan Tianchi Volcano Observatory, are plotted in Figure 8. Additionally, the leveling measurement by Hu et al. [36] and the surface deformation results estimated by Tang P. et al. [35], He et al. [32], and Kim et al. [34] are plotted in Figure 8.

Figure 8 shows that the maximal magnitude of the earthquakes linearly increases from M_L 2.7 in 2002 to M_L 4.4 in 2004 but slightly decreases to M_L 4.0 in 2005, while the number reached highest in 2003. Meanwhile, the CSD shows a good agreement with the maximal magnitude of the earthquakes from 2002 to 2004, but the agreement is not good from 2004 to 2005.

From 2007 to 2010, the retrieved CSD is consistent with the change of earthquake magnitude and frequency, but the deformation range is quite different. E.g., the CSDs by Tang P. et al. [35], He et al. [32], and Kim et al. [34] are consistent with maximal magnitude of the earthquakes. Particularly, the variation of the maximal magnitude of the earthquakes is clearly lower from 2008 to 2009 than that from 2009 to 2010, and the CSD shows a similar behavior as such.

Combined with the findings mentioned above, the time series CSD is strongly correlated with the seismic activities in the Changbaishan volcano, where the annual maximal earthquake magnitude likely dominates the surface deformation in the corresponding year.

4.3. Re-Evaluating Deep Pressure Source of Changbaishan Volcano

In the past 20 years, the location and geometric characteristics of the deep pressure source are always the hotspots in studying the Changbaishan volcano [8,9].

4.3.1. Model Selection

Hu et al. [58] analyzed the possible shape of deep pressure source by comparing the effects of the single and double Mogi source models, prolate spheroid, and Okada rectangular dislocation model on the surface deformation in the Changbaishan volcano. Hu et al. [36] inversed the magma chamber with single Mogi source model using the GPS and leveling data from 2002 to 2005, and found that the magma chamber was shrinking. Chen et al. [30] studied the magma chamber of the Changbaishan volcano based on the single Mogi and double Mogi source models with the leveling measured and GPS data. Zhu et al. [59] inverted the deep pressure source with GPS and leveling measured data from 2002 to 2003 using three mutually perpendicular expansion point source models. He et al. [32] assessed the magma volume change rate using single Mogi source model with the InSAR surface deformation from 2006 to 2011. Kim et al. [34] simulated a deep pressure source combined with the single Mogi source and Okada dislocation models. Huang et al. [9] used the finite element equivalent physical strength method to invert the deep pressure source of the Changbaishan volcano, which indicated that the surface deformation caused by an ellipsoidal magma chamber fits well with the observed GPS and leveling measured data. By defining a spheroid source, a tabular reservoir, and a NW-SE fault, Trasatti et al. [42] indicates that the Changbaishan volcano is still affected by deep magmatic recharge from 2018 to 2020, and a NW-SE trending fault is reactivated. These results about the deep pressure source of the Changbaishan volcano are summarized in

Table 2. Previous inversion results of deep pressure source of the Changbaishan volcano.

	Model	Depth (km)	ΔV$(\times {10}^6 {\mathbf{m}}^3)$	Strike (°)	Plunge (°)
Hu et al. [58]	Double Mogi source	3.2791 (2002–2003)	0.004	-	-
		4.1968	0.015	-	-
Hu et al. [36]	Single Mogi source	6.854 (2002~2003)	16.980	-	-
		7.653 (2003~2004)	10.620	-	-
		3.1515 (2004~2005)	3.390	-	-
Zhu et al. [59]	Okada dislocation	9.2 (2002–2003)	5.220	337.7	14.0
		9.2	7.0	34.4	82.2
		9.2	6.598	302.8	78.5
Chen et al. [30]	Double Mogi source	7.9 (2002–2003)	6.3	-	-
		5.5	1.6	-	-
He et al. [32]	Single Mogi source	5.3 (2006–2009)	2.475 *	-	-
		6.9 (2009–2011)	2.093 *	-	-
Kim et al. [34]	Single Mogi source	11.465 (2007.8–2007.12)	3.756	-	-
		11.094 (2008.3–2008.5)	1.579	-	-
		10.013 (2010.0–2010.4)	−1.557	-	-
Huang et al. [9]	Finite Element	6.86 (2002–2003)	-	77.5	44.3
			-	324.5	69.1
			-	217.9	53.1
Trasatti et al. [42]	Prolate spheroid	7.7	7.3 (±2.4)	140	0
	Tabular reservoir (2018–2020)	14.0	−14 (±4)	50	0
	Single Mogi source (2002–2005)	4.40	4.6 (±1.53)	-	-
	Single Mogi source (2009–2011)	6.70	−1.4 (±0.2)	-	-

Note: The source volume retains two decimal places, and the data marked with * are converted from the data in the original paper.

.

According to the previous results, we thought that the shape of the deep pressure source of the Changbaishan volcano is likely irregular. Therefore, to evaluate the deep pressure source of the Changbaishan volcano, we have selected four idealized deep pressure source models with different geometric shapes, including the single Mogi [60], penny-shaped sill-like [61], finite spherical [62], and prolate spheroid [63] sources. Moreover, considering that deep pressure sources may distribute in different depths and there probably exists more than one layer of magmatic sources [7,8,32,64], the double Mogi source model is also introduced in this study.

4.3.2. Inversion of Volcanic Deep Pressure Source

To evaluate the deep pressure source of the Changbaishan volcano, the Bayesian algorithm framework is adopted [65]. First, the CSD is simulated using the aforementioned five deep pressure source models, and the likelihood function estimation is obtained by comparing the simulated CSD and the inverted CSD in this study when considering the errors of the Sentinel-1B InSAR data. Then, based on the prior information and likelihood function estimation, the Markov chain Monte Carlo algorithm and Metropolis-Hastings algorithm are used to iteratively sample the posterior probability density function of deep pressure source model parameters. Finally, through finding the maximum posteriori probability solution of the posteriori probability density function, a set of optimal deep pressure source model parameters (e.g., the horizontal position, depth, and geometric structure) can be calculated.

Here, the Bayesian algorithm is operated with the MATLAB^® software package, GBIS (Geodetic Bayesian Inversion Software) V1.1 [65]. Since the variation of the maximal magnitude and the number of annual earthquakes is relatively slight from 2007 to 2020 (Figure 8), the prior information of single Mogi and double Mogi source models refers to the relevant researches from 2007 to 2016 [32,34]. The prior information regarding the penny-shaped sill-like, finite spherical, and prolate spheroid source models are referred to in previous studies [9,42,59]. The prior information of the five deep pressure sources is present in

Table 3. The prior information for the single Mogi, double Mogi, penny-shaped sill-like, finite spherical, and prolate spheroid sources.

Model	Parameter	Start	Step	Lower	Upper
Single Mogi and double Mogi	X (m)	0	50	−15,000	15,000
	Y (m)	0	50	−15,000	15,000
	Depth (m)	3000	100	100	10,000
	DV (m3)	1 × 105	1 × 105	−1 × 107	1 × 109
Penny-shaped sill-like	X (m)	0	100	−20,000	20,000
	Y (m)	0	100	−20,000	20,000
	Depth (m)	2000	50	100	20,000
	Radius (m)	1000	100	100	2000
	ΔP/μ	−1 × 10−3	1 × 10−4	−1 × 103	1 × 103
Finite spherical	X (m)	−1000	100	−20,000	20,000
	Y (m)	1000	1000	−20,000	20,000
	Depth (m)	3000	100	1000	15,000
	Radius (m)	1000	50	100	5000
	ΔP/μ	−1 × 10−3	1 × 10−4	−1 × 106	1 × 105
Prolate spheroid	X (m)	0	100	−20,000	20,000
	Y (m)	0	100	−20,000	20,000
	Depth (m)	5000	50	1000	20,000
	Major semi axis (m)	500	50	1	20,000
	Aspect ratio	0.5	0.01	0.01	0.99
	ΔP/μ	−1 × 103	1 × 10−4	−1 × 103	1 × 103
	Strike (°)	270	1	1	359
	Plunge (°)	−30	1	−89.9	89.9

Note: The coordinates X (km) and Y (km) are centered in the geo-reference point: [128.084°, 42.051°]; Depth (km) is depth of deep pressure source (positive downward); DV (×106 m³) is the change of the deep pressure source volume; Radius (m) is the radius of the sphere; Major semi axis (m) is the length of major semi-axis; Aspect ratio is dimensionless aspect ratio between semi-axes (minor/major); ΔP/μ is the dimensionless excess pressure; Strike (°) is the strike angle of major semi-axis with respect to North (0°/360° = N; 90° = E; 180° = S; 270° = W); Plunge (°) is the inclination angle of major semi-axis with respect to horizontal (0° = horizontal; 90° = vertical).

.

To improve the efficiency of Bayesian inversion, the simulated annealing algorithm is used to obtain the “Start” parameters of each deep pressure source model in the initial stage of inversion. Based on the above priori parameters and the Bayesian algorithm, the optimal posteriori probability solution, mean, median, and 2.5 percent and 97.5 percent of posterior probability density functions of location, depth, volume, and other parameters derived with the aforementioned five deep pressure source models are obtained (

Table 4. The ‘Optimal’, ‘Mean’, ‘Median’ values, and 97.5% confidence of inverted deep pressure source parameters.

Model	Parameter	Optimal	Mean	Median	2.5%	97.5%
Single Mogi	X (km)	−5.315	−7.126	−6.910	−10.809	−4.603
	Y (km)	3.172	5.946	5.899	2.126	9.967
	Depth (km)	6.2	6.2	6.4	3.8	7.7
	$\mathrm{DV} (\times$106 m3)	3.2	5.7	5.3	2.5	10.6
Double Mogi1	X (km)	−3.831	−3.817	−3.847	−5.109	−2.852
	Y (km)	1.447	1.913	1.689	0.745	4.330
	Depth (km)	4. 7	4.7	4.8	3.0	6.3
	$\mathrm{DV} (\times$106 m3)	1.8	2.0	2.0	0.9	3.8
Double Mogi2	X (km)	−7.088	−7.192	−7.124	−7.412	−7.075
	Y (km)	−0.869	−0.909	−0.869	−1.009	−0.334
	Depth (km)	0.012	0.106	0.042	−0.004	0.334
	$\mathrm{DV} (\times$106 m3)	0.0006	0.0171	0.0017	0.0007	0.0658
Penny-shaped sill-like	X (km)	−0.741	−0.931	−0.938	−1.180	−0.677
	Y (km)	−5.908	−5.665	−5.655	−5.968	−5.358
	Depth (km)	1.4	1.9	1.9	1.3	2.1
	Radius (m)	236.2	341.0	332.1	235.3	427.6
	ΔP/μ	−0.01	−0.00	−0.00	−0.01	−0.00
Finite spherical	X (km)	3.304	3.380	3.372	2.971	3.834
	Y (km)	2.428	2.491	2.476	1.966	3.098
	Depth (km)	1.0	1.0	1.0	1.0	1.2
	Radius (m)	149.5	227.3	170.2	102.0	686.1
	ΔP/μ	−0.56	−0.60	−0.39	−1.86	−0.01
Prolate spheroid	X (km)	0.618	0.687	0.685	0.382	1.041
	Y (km)	−0.708	−0.741	−0.719	−1.195	−0.413
	Depth (km)	1.0	1.1	1.0	1.0	1.3
	Aspect ratio	0.04	0.04	0.02	0.01	0.13
	ΔP/μ	−0.04	−0.18	−0.11	−0.67	−0.00
	Major semi axis (m)	2378.4	2355.7	2357.1	1879.0	2838.6
	Strike (°)	291.18	291.57	291.64	281.12	312.12
	Plunge (°)	53.18	53.19	53.51	46.29	57.03

Geo-reference point: [128.084°, 42.051°].

).

4.3.3. Inversion Results

Table 4 shows that, among the five volcanic deep pressure source models, the estimated “Depth” with the single Mogi source model is 6.2 km, which is closer to the results of Hu et al. [36] (6.854 km in 2002 to 2003), He et al. [32] (6.9 km in 2009–2011), Huang et al. [9] (6.86 km in 2002 to 2003), and Trasatti et al. [42] (6.7 km in 2009 to 2011, 7.7 km in 2018 to 2020) in Table 2. In addition, the location (X, Y) of the single Mogi source is located in the northwestern slope of the Changbaishan volcano, which is consistent with the viewpoint of Yang et al. [64], who proposed that the deep magma chamber of the Changbaishan volcano is not directly below the Changbaishan volcano, as well as with the conclusion of He et al. [32], who given that the single Mogi source is located in the west of the lower part of the volcano. Therefore, the inversion result with single Mogi source model should be rational.

Then, the posterior probability density functions of the deep pressure source parameters can be obtained through 1 × 10⁷ iterations of the single Mogi source model, as is shown in Figure 9. The scatter diagrams in the first three rows show the joint distributions between the parameter pairs, and the bottom rows show the histograms of marginal distributions of each parameter. The warm color represents the parameters with high frequency, corresponding to the horizontal and vertical axis, while the cold color is the opposite. Figure 9 intuitively reflects a good correlation between the parameter pairs, using the single Mogi source model, which again validates that the single Mogi source parameters are reasonable. Here, the optimal estimated depth of the single Mogi source is 6.2 km, and the volume increase is 3.2 × 10⁶ m³. Compared with the prolate spheroid source volume variation of Trasatti et al. [42] in 2018 to 2020, our single Mogi source results (optimal parameters) are consistent with it in order of magnitude. The difference between the two is 4.1 (±2.4) × 10⁶ m³. Considering that we have monitored for a longer period (2016–2020), and the variation of surface deformation will make our results different from it, which is reasonable.

Therefore, the complex local surface deformation in the Changbaishan volcano can be partly explained by the expansion of the single Mogi source distributed in the northwest of Tianchi caldera. In the time series, it shows a strong correlation with the maximal magnitude and number of annual earthquakes. Additionally, in space, it shows that the locally distribution of the CSDs is bounded by faults, and the largest uplift mainly occurs in the western slope of the Changbaishan volcano. In recent years, earthquakes frequently occurred in the Changbaishan volcano. On 5 March 2021, an M_L 3.1 earthquake has occurred in Fusong, China (41.81°N, 127.87°E, depth 25 km, https://www.cea.gov.cn/cea/xwzx/zqsd/5576146/index.html, accessed on 6 March 2021), whose epicenter is located just along fault f2. This supports the finding that fault f2 plays an important role in the time series CSDs in the Changbaishan volcano. Therefore, it is inferred that the tectonic activity in the Changbaishan volcano is active, and the faults related to the deep pressure source are also in the active state.

5. Conclusions

In this study, the surface deformation of the Changbaishan volcano in the LOS direction is retrieved using the SBAS-InSAR method with 75 Sentinel-1B SAR data, from 2016 to 2020. The deep pressure source is evaluated with five models, based on the Bayesian algorithm via the calculated CSD. The main findings are presented as follows:

(1)

The MDV in the Changbaishan volcano is uplifted as a whole. The uplifted extent is relatively high in the east part, relatively low in the middle region, located in the eastern and southern slopes, but highest in the western and northern slopes.

(2)

There exists a seasonal fluctuation in the time series CSD of the Changbaishan volcano. It shows subsidence mainly from October to February or March of the following year and uplift from April to August or September.

(3)

The CSD shows a strong correlation with the maximal magnitude and the number of annual earthquakes. After comparison, combined with the previous results, the CSD in the Changbaishan volcano is likely dominated by the maximal magnitude of the earthquakes in the corresponding year.

(4)

There is a strong correlation between surface deformation and the faults in Tianchi Volcano. Fault f2 plays a controlling role on the surface deformation in the Changbaishan volcano, followed by faults f6 and f1.

(5)

After comparison, the single Mogi source model is appropriate for evaluating the deep pressure source in the Changbaishan volcano, constrained by the calculated surface deformation in this study. The results show that the optimal estimated depth of the magma chamber is 6.2 km, and the volume increase is 3.2 × 10⁶ m³. Additionally, it is inferred that the tectonic activity in the Changbaishan volcano is in an active state.