Changbaishan intraplate volcano (2,774 m a.s.l., also known as Paektu, Baekdu, Tianchi; Figure 1A) is located at the China/North Korea boundary above the 500–600 km deep stagnant portion of the Pacific slab (Zhao, 2021). The volcano is affected by a summit caldera depression with a diameter of about 5 km filled by a lake, and its activity developed in three main stages (Zhang et al., 2018): a shield-forming stage (5.02–1.05 Ma), a cone-construction stage (1.37–0.01 Ma and mostly<1 Ma), and the more recent, caldera-forming stage (0.01 Ma to 1903 CE). This last stage produced trachytic-comenditic pumices during Plinian eruptions. The 4.77 Ma to 1.88 Ma old central volcanoes Wangtian’e and Namphothe are located about 40 km southwest and southeast of Changbaishan, respectively. The basement of Changbaishan is constituted by the Gaima basalt plateau and the underlying crystalline and metamorphic rocks. The Gaima plateau and the Chagbaishan flanks are characterized by numerous monogenetic scoria cones roughly aligned along a NW-SE direction (Zhang et al., 2018; Figures 1A,B). This is also the strike of the active Baishan-Jince fault, the main structural discontinuity affecting the basement and Changbaishan edifice, and along which crustal earthquakes and monogenetic cones align (Liu et al., 2004; Figures 1A,B). The Changbaishan morphology consists of a slightly NW-SE elongated cone with a sub-circular caldera mainly formed during the M 7.2, “Millennium” eruption (ME) in 946 ± 20 CE (Figure 1) (Zhang et al., 2018). A NW-SE striking, SE dipping fault affects the western slope of the cone (Liu et al., 2004; Figure 1). Except for the earthquakes recorded during the 2002–2006 unrest, which concentrated in the Changbaishan caldera, the 1977–2018 events roughly align along this fault (Figure 1A). The above summarized structural, geophysical, and morphological data coherently indicate that the volcano formed on a NW-SE weakness zone possibly controlled by the Baishan-Jince regional fault system, which extends from North Korea to China. Changbaishan represents one of the most dangerous Eastern Asia volcanoes because of its ‘Millennium’ eruption (ME), which emitted about 96 ± 19 km³ of pyroclastics and whose gas emissions had a global impact on climate (Horn and Schmincke, 2000; Iacovino et al., 2016; Yi et al., 2021). The 2002–2006 unrest characterized by deformations, earthquakes including low-frequency events, and changes in volcanic gas geochemistry (Xu et al., 2012). CO₂-rich degassing areas and hot water are widespread in the Changbaishan area with CO₂ discharge average values of 1.8 × 10⁵ t/yr (Zhang et al., 2015). The ³He/⁴He of the emitted gas (Ra = 6.6 in 2006) and the B isotopic values suggest that the Changbaishan fluids are released by the asthenospheric mantle (Xu et al., 2012) metasomatized by recycled carbonates and organic metasediments of the Pacific stagnant slab (Zhang et al., 2015).

Thermobarometry petrological data indicate that the plumbing system consists of sill-like reservoirs located at different depths (Lee et al., 2021; Yi et al., 2021). Seismological, resistivity and gravity data (Choi et al., 2013; Qiu et al., 2014; Hammond et al., 2020) evidence the possible occurrence of melts between 4 and 8 km depth while the hydrothermal system is located at about 3 km depth (Liu et al., 2021). Petrological studies on the evolution of the Changbaishan magmatism based on clinopyroxene thermobarometers (Lee et al., 2021) indicate a large, sill-like reservoir located between 10 and 13 km depth, while fluid inclusion studies (Andreeva et al., 2018) confirm this deep reservoir and evidence, according to the previously reported geophysical data, two smaller and shallower reservoirs at 7–8 km depth and at 3.5–5 km depth. During the 2002–2006 unrest, the hypocentral depth of earthquakes ranged between 12 km and the surface with most events concentrated in the upper 4–5 km of the crust below the caldera (Liu et al., 2021). A strong reduction in the number of earthquakes occurred in the range 5–7 km and 9–12 km suggesting partly molten zones within the present-day plumbing system.

GNSS, InSAR and leveling data do not record significant deformations in the 1992–98 period, whereas a 6.5 cm maximum uplift has been detected during the unrest from 2002 to 2006 (Ji et al., 2021). This uplift was followed by a subsidence between 2008 and 2011. By using GPS and levelling data and by applying a Mogi-like source, Xu et al. (2012) suggest that the depth of the 2002–2006 inflation was at 6 km, whereas Zhu et al. (2008) propose a depth of 9.2 km. However, Ji et al. (2021) remark that the low spatial resolutions and geographic distribution of the geodetic sites do not allow to better constrain the source of the 2002–2006 uplift. This issue was partly overcome by InSAR data, whose larger spatial coverage constrains a source located at about 9 km depth for the 2002–2006 unrest (Ji et al., 2013). Chen et al. (2008) analyze GPS, levelling, ENVISAT SAR and JERS-1 SAR data collected in 2002 and 2003 and conclude that the deformation field is consistent with two Mogi-type pressure sources: a deeper source at 7.9 km depth located below the Changbaishan caldera and a shallower one at 5.5 km depth located about 40 km southeast of the caldera in North Korea below the Namphothe volcano (Figure 1). These two sources are connected by the Baishan-Jince regional fault system. The post-unrest subsidence and the more recent dynamics of Chagbaishan are poorly known, lacking specific studies and available data. The only available data concern the seismicity (Pan et al., 2021); these data show that after 2002–2005, when an average of 500 events/yr with M_Lmax = 4.4 was detected, the seismicity sharply decreased to 30 events/year with M_Lmax = 2.8 between 2010 and 2018.

The Interferometric Synthetic Aperture Radar (InSAR) technique allows the estimation of ground motion along the sensor Line-Of-Sight (LOS) direction. Traditional differential InSAR (DInSAR) consists of a differential interferogram between two co-registered SAR scenes acquired at different times over the same area, also considering a digital elevation model (DEM). The main limitations are due to geometrical-temporal decorrelation and atmospheric artifacts (e.g., Peltzer and Rosen, 1995). Targeting at more precise deformation monitoring, several multi-temporal InSAR (MTI) techniques have been proposed in the last decades, which exploit the redundancy offered by hundreds of image pairs to reduce the above-mentioned limitations (e.g., Ferretti et al., 2001; Berardino et al., 2002). The outcomes of MTI consist of a mean ground velocity map and the relative displacement time-series temporally referred to the first acquisition date of the image stack, and spatially referred to a stable reference point or area within the radar coverage.

One of the most adopted MTI algorithms is the Persistent Scatterers (PS) InSAR approach that is a proven differential interferometric technique which involves processing of multitemporal SAR data to identify persistently scattering ground features and their motion rates with millimeter precision (Ferretti et al., 2001). PS concept is based on finding permanent scatterers showing phase stability over a long period of time and aims to measure accurate surface movements. The advantage of such a technique is that it overcomes geometrical and temporal decorrelation issues in urban and rocky areas by using a large amount of radar images. The PS analysis starts with a linearly independent set of interferometric combinations with respect to the same master scene. After removal of the topographic phase using the simulated phase relative to a DEM, the geocoded pixels are further investigated based on their phase stability. Afterwards, PS candidates with stable phases are selected with a given threshold. To account for residual orbit errors, phase ramps can be calculated and removed from the interferograms using a set of ground control points. Finally, atmospheric effects can affect the remaining phase of the interferograms, so a double spatio-temporal filtering of the resulting time series is performed to separate the deformation from the atmospheric signals.

Due to the problem of heavy vegetation coverage, rugged topography, long time snow coverage (from October to May) and cloud coverage of the study area, the C-band (ENVISAT and Sentinel-1) SAR images collected by European Space Agency (ESA) suffer from extremely low coherence. On the contrary, the L-band ALOS-2 data allows to overcome such issues providing sufficient coherence and spatial density of PS points. Therefore, we process 19 ALOS-2 images (path 131, frame 830) from ascending orbit in the StripMap acquisition mode, provided by the Japanese Aerospatial Agency (JAXA). The incidence angle is 40.58° and the temporal spanning is between 20th November 2018 and 20th October 2020.The multi-temporal algorithm named Persistent Scatterers was adopted with the ENVI SARscape software (L3HarrisTM, https://www.l3harris.com/). Also, the April 9, 2019 acquisition is selected as master image (Supplementary Figure S1). The Shuttle Radar Topography Mission (SRTM) (Farr et al., 2007) with 1 arc-second (∼30 m) pixel size DEM is used to remove the topographic phase contribution. Due to cloud coverage and rugged topography of Changbaishan area, the tropospheric phase error, which includes vertical stratification and turbulence, could affect the accuracy of deformation estimation (Parker et al., 2015). Here, the vertical stratification is removed using the SRTM DEM. For the removal of tropospheric turbulence, both a low-pass filtering of 1,200 m and a high-pass one of 365 days are performed. Finally, a ground pixel of 30 m is considered during the geocoding step. Snow coverage and dense vegetation are both important factors of coherence loss in multi-temporal InSAR. However, L-band SAR images, such as the ALOS-2 ones, have the characteristic to penetrate sparse vegetation and snow coverage. There are only 4 months (from June to September) without snow coverage on the Changbaishan volcano summit. Also, by comparing the coherence with the land cover (Ji et al., 2010), dense vegetation (e.g., forest, mixture of more forest and less soil) is the main reason for coherence loss in Changbaishan area. A coherence threshold of 0.8 is chosen to select the final set of PS, and areas around the caldera provide a sufficient amount of PS points.

Initially, InSAR data have been downsampled considering a varying length increasing from 250 m close to the crater, to 500 m and up to 1,500 m in the far-field. The inversion of geodetic data is performed by a Bayesian global optimization tool. The Bayesian inference is one of the inversion options available in VSM—Volcano and Seismic source Modelling (Trasatti, 2019). VSM implements the most common deformation sources that may be combined, and the most common geodetic data-types (SAR, GNSS, leveling, etc.). The Bayesian optimizer samples the parameters space using a Monte Carlo setting and combines the a-priori information with the likelihood computations to obtain posterior probability density functions of the parameters (Mosegaard and Tarantola, 1995). The Changbaishan active sources are represented by the above defined plumbing system and active fault. The medium is assumed as homogeneous and isotropic with shear modulus μ = 30 GPa and Poisson ratio equal 0.25. Inelastic rheology is likely to take place at volcanoes; however due to the limited time window covered and the seismicity located in the upper crust below the volcano, we consider mean ground velocities and linear elasticity of the medium. After several attempts, we choose an inflating spheroidal source (Yang et al., 1988) instead of an isotropic point-source (Mogi, 1958) to reproduce the NW-SE elongated high positive LOS velocities. The assumed source is a NW-SE elongated horizontal spheroid with fixed strike N140°E and axes proportion 1:3. A second volumetric source characterized by negative volume variation is fixed at 14 km b.s.L. and it is represented by a square tabular reservoir undergoing negative perpendicular displacement (Okada, 1985). This deep reservoir is centered below the caldera. Then, a NW-SE dip-slip fault is considered (Okada, 1985), whose characteristics such as strike, dip, position, and extension are based on previous studies (Figure 1) and on the InSAR velocity pattern depicted in Figure 2. We are aware that different settings of the active sources may be equally consistent with data. Our model configuration of the plumbing system is compatible with data and, at the same time, limits the number of unknowns.

The mean parameters resulting from the inversion of the InSAR data are reported in Supplementary Table S1, while the synthetic velocities and their comparison with data are reported in Figure 3. The 1D and 2D posterior probability density functions of the inverted parameters are reported in the Supplementary Figure S2, along with the samples collected during the random walks (Supplementary Figure S3). Comparisons with subsampled data are reported in Supplementary Figure S4.

The spheroid lays SE of the caldera, at 8.7 km below the reference height here fixed at 1,000 m, therefore we consider it to be located at about 7.7 km b.s.L. The reference height is usually chosen as the average elevation of the used data in case of remarkable topography (Bonaccorso et al., 2005). The volume variation estimated for the inflating spheroid is 7.3 ± 2.4 106 m³/yr, while the retrieved finite volume amounts to about 300 km³. The volume variation of the deeper tabular reservoir (laying 14 km b.s.L.) is -14 ± 4 106 m³/yr. The amount of contraction of the deep reservoir in terms of its decreasing annual volume is almost twice the inflation of the shallower inflating source.

Figure 3 reports the results of the modeling and comparisons with data. In particular, the NW-SE trending main uplift area is reproduced by the spheroidal source. Also, the presence of the active dip-slip fault inhibits inflation in the western sector of Changbaishan volcano, as shown in Figure 3B. The far-field displacements are also reproduced on average. Figure 3C shows residuals, i.e., the difference between data and modelled LOS displacements. Pixels in white are within the uncertainty associated with InSAR data, here fixed at 4 mm/yr according to, e.g., Casu et al. (2006). Figures 3D,E shows the comparison between data and model along the NW-SE and SW-NE profiles depicted in Figure 3A. The ALOS-2 velocities are quite scattered, and the model reproduces them on average. The negative to positive velocities along the NW-SE profile are reproduced by the combined action of the shallow inflating source and the deep tabular deflating source. The peaks on the SE are due to residual tropospheric turbulence and cannot be reproduced by the analytical models. On the other hand, the SW-NE profile clearly highlights a jump in the velocities while crossing the NW-SE trending fault (grey vertical line in Figure 3E). The negative velocities on the western side of the volcano, west of the fault, and the short-wavelength positive velocities moving towards the caldera are simulated by the dip-slip mechanism combined with the spheroid inflation. We exclude that the deformation pattern depicted in Figure 3E (SW-NE profile) could reflect magma related dynamics including cooling and/or degassing because there is no evidence of magma reservoirs in the western sector of the volcano (Ji et al., 2010, 2013; Choi et al., 2013; Hammond et al., 2020). The deformation field is consistent with movements along an outcropping fault. This fault may move in response to gravity-driven processes on the western flank of the volcano. This is supported by the steep topography (Figure 3E) and by the transition from hard rock (lava flows and domes) to partly loose material (pyroclastics). Finally, data close to the caldera rim show local minima/maxima and may be affected by residual geometrical artifacts in the InSAR analysis, due to the steep slopes and/or localized gravity processes.

This study is aimed at comparing the features of the plumbing system active during the 2000s with those retrieved for the present unrest. We consider the literature leveling data related to the 2002–2005 unrest and successive 2009–2011 subsidence (data redrawn from Xu et al., 2012). The 2002–2005 unrest shows a maximum of 6.5 cm uplift close to the caldera along the North route (levelling routes shown in Figure 2). Since 2006, the ground uplifted at a slower rate until 2009, when a subsidence phase started. For the subsidence phase we consider elevation changes from North and West leveling routes, amounting to -8 mm maximum close to the caldera.

The Bayesian inversions are performed using the above mentioned VSM tool. To limit the number of unknown parameters due to the few available benchmarks, we consider only part of the deformation sources investigated above. The 2002–2005 leveling data consist of the North levelling route alone; therefore we consider only a simplified isotropic point-source (Mogi, 1958). The horizontal position of the isotropic point-source is kept equal to that obtained for the spheroidal source, while its depth and volume variation may vary. Since the 2009–2011 data account for the western sector also, we include the NW-SE fault with parameters fixed as for the inversion of the InSAR data (Supplementary Table S1). Results of the leveling data inversions are summarized in Supplementary Table S2 and Figure 4.

The uplift during the 2002–2005 unrest is best modelled by a point-source placed at 4.4 km b.s.L., with a volume variation of 14 × 106 m³ accumulated in 3 years, i.e., an annual rate of 4.7 × 106 m³/yr. The source is shallower than that retrieved by ALOS-2 data and its volume variation annual rate is less than that obtained with ALOS-2 data. For the 2009–2011 subsidence, we retrieve a point-source at 6.7 km b.s.L., slightly shallower than that retrieved by InSAR. The annual contraction amounts to -1.4 × 106 m³/yr, a small fraction of that retrieved for the spheroid depicted by InSAR. The NW-SE fault undergoes a dip-slip of about 1 mm/yr, a value much smaller than that retrieved by InSAR. These values may be considered indicative because of the limited resolving power of the few leveling measurements.

A sketch of the retrieved sources at Changbaishan is reported in Figure 5. The results of our analysis of the 2018–2020 ALOS-2 data show that Changbaishan is affected by a general uplift of its southeastern flank with values up to 20 mm/yr and by a significant subsidence, down to –20 mm/yr of its southwestern sector at altitudes less than about 2000 m. The uplift zone is NW-SE elongated and mainly extends in the North Korea sector of the edifice, where the maximum deformation is recorded. A similar deformation field has been detected by Chen et al. (2008) in 2002 and 2003 from the analysis of ENVISAT SAR and JERS-1 SAR data. The maximum elongation of the 2018–2020 uplifting area detected by our analysis follows the Baishan-Jianchi fault regional fault system. Also, the results of our modeling indicate that the NW-SE striking fault affecting the western flank of the volcano separates the zone of uplift from that of subsidence. The occurrence of past earthquakes along this discontinuity (Figure 1) and our deformation data reveal that this fault is active and moves by creep with normal mechanism and slip-rate of 0.019 m/yr during 2018–2020. This slip-rate corresponds to a geodetic scalar moment in the order of 10¹⁵ Nm/yr, equivalent to the energy released by a Mw = 3.9 earthquake (Kanamori, 1977). The creep movement of faults during the inflation/deflation phases is not peculiar to Changbaishan but has been detected at other volcanoes, e.g. at Campi Flegrei (Vilardo et al., 2010), Etna (Bonforte et al., 2007) and Ischia (Trasatti et al., 2019), Italy.

At Changbaishan, the evidence of the 2018–2020 NW-SE elongated nearly-inflation and normal fault-like pattern suggests that the deformation field is possibly controlled by the NW-SE striking Baishan-Jianchi fault system (Figure 1). According to Zhang et al. (2018), most of the monogenetic cones concentrate along this fault system, suggesting that it represents a preferential path for magma ascent (Figure 1B). Our results of the modeling of the recorded uplift suggest that a tabular, deflating deep source at 14 km b.s.L. decreases its volume at a rate of about −14,106 m³/yr, while the NW-SE elongated, shallower source at 7–8 km depth inflates at a rate of 7.3 106 m³/yr (Supplementary Table S1). The geometry and depth of such sources are consistent with the results from petrological data (Andreeva et al., 2018; Yi et al., 2021; Lee et al., 2021) and gravity modeling (Choi et al., 2013), which show that the Changbaishan plumbing system consists of magma reservoirs with tabular or sill-like geometry located at different depths. Andreeva et al. (2018), Lee et al. (2021) and Choi et al. (2013) suggest the presence of a deep reservoir at about 10–15 km depth in which basaltic melts easily accumulate and upraise to a shallower, from 8 km to about 5.5 km deep reservoir where magma may evolve to form basaltic trachyandesites and trachytes. The low number of earthquakes below 10 km depth and between 3 and 7 km depth (Liu et al., 2021) supports the hypothesis of the occurrence of such magma reservoirs at these depths (Figure 6). Receiver function (Hammond et al., 2020) and magnetotelluric data (Qiu et al., 2014) also reveal a reduction of seismic velocities and a low resistivity zone between 4 and 8 km depth, suggesting the occurrence of melt. Considering these independent estimates regarding the depth and geometry of the Changbaishan magma reservoirs, we conclude that the tabular and deflating source at 14 km depth evidenced by our modeling of the ALOS-2 data is compatible with the 10–15 km deep basaltic reservoir. Also, the shallower, inflating spheroidal reservoir at 7.7 km b.s.L. may correspond to the 4–8 km deep basaltic trachyandesitic to trachytic magma chamber. Then, the NW-SE elongation of the spheroid reflects the tectonic setting of the region since the magma storage takes place in optimal stress orientation regime. Therefore, the 2018–2020 deformation field records a transfer of magma from a deeper to a shallower source, a mechanism observed at other volcanoes, e.g. at Grímsvötn, Iceland (Reverso et al., 2014).

Here, we consider the inferred volumetric changes deduced for the two reservoirs and use their difference to determine the bulk modulus of the magma in the reservoir. The relation linking the volume change of a spherical magma chamber ∆V _c and the magma flow ∆V _m into/out is (Johnson, 1992; Johnson et al., 2000; Hautmann et al., 2017): Δ V m     =   2 3   ( 4 μ 3 K     ) Δ V c   with μ as the crustal shear modulus and K as the bulk modulus of the magma. By using |∆V _c | = 14 × 106 m³/yr and ∆V _m = 7.3 × 106 m³/yr (Supplementary Table S1) and adopting an average crustal shear modulus of 30 GPa, we determine a magma bulk modulus of 51 GPa, a value consistent with that of basaltic melts (K = 50–60 GPa between 1,000 and 1,200°C; Schilling et al., 2003). This should be considered in terms of order of magnitude estimate because of the spherical geometry assumption. Considering this result, the volume deficit observed in the transfer of magma from the deflating and deeper reservoir to the shallower one cannot be attributed to magma compressibility (=1/K), i.e. magma plus gas. Therefore, we suggest, according to the crystal-mush model proposed by Yi et al. (2021) for the plumbing system of Changbaishan, that about half of the gas-poor melt volume exiting from the 14 km deep reservoir intrudes a porous crystal-mush, while the remaining half enters the 7.7 km deep sill-like reservoir. The conduit separating the two reservoirs has a vertical length of about 6.3 km. By ignoring horizontal extensional stresses and assuming an overpressure P _f = 10 MPa (Manga and Brodsky, 2006), the width a of a dike-like conduit may be calculated by the relation (Hartley and Maclennan, 2018): a     =   h P f /   μ by selecting μ = 30 GPa and h = 6.3 km, we obtain a ∼ 2 m. The vertically directed upward force on the magma column between the two reservoirs is given by ∆ρgh, with ∆ρ the density difference between magma and host rocks, g the gravity acceleration and h the difference between the depths of the two reservoirs. Assuming ∆ρ between 100 and 200 kg/m³, the vertically directed upward force on the column is about 13–27 MPa, a value larger than the 1–10 MPa pressure required for basaltic reservoirs to open a crack so allowing magma uprising (Manga and Brodsky, 2006). We underline that the volume variation of the conduit amounts to 2 × 104 m³, that is 2-3 orders of magnitude less than the volume variation of the two reservoirs, and therefore it is not constrainable by geodetic data.

Our results from the inversion of levelling data (Xu et al., 2012) relative to the 2002–2005 uplift during the unrest phase and the 2009–2011 subsidence consist in:

1) a point-source inflating between 2002 and 2005 with a volume variation rate of 4.6 × 106 m³/yr and depth 4.4 km b.s.l. As previously reported, the source is shallower than that retrieved by ALOS-2 data and the annual rate of volume variations is also less than that recorded for the 2018–2020 period. This is due to the shallower depth retrieved for the 2002–2005 source. Both the 2002–2005 point-source and the 2018–2020 spheroid are coaxial.

The above summarized values must be considered in terms of their order of magnitude, because, according to Ji et al. (2021), the limited distribution of the geodetic sites on the northern and western sides of the volcano does not allow us to constrain the whole deformation pattern of Changbaishan. However, some inferences on the processes occurred during the uplift and the following subsidence phase may be suggested. The increase in CO₂ concentration and values of the He isotopic ratio (Xu et al., 2012) in the emitted gases during the uplift indicate the involvement of a pressurized degassing magmatic source possibly located, according to our modeling results, at 4.4 km depth (Figure 6). This conclusion is supported by the lack of earthquakes at this depth during 2002–2006 (Liu et al., 2021). While the 2002–2006 uplift mainly recorded the arrival of gas to the top of the reservoir, the successive 2009–2011 subsidence, whose source is about 2 km deeper than of the uplift, can be associated to gas release and crystallization processes in the middle/lower portions of the reservoir. This mechanism has been proposed to explain the deflation following previous uplift phases associated with unrest episodes at other volcanoes, e.g., at Okmok volcano, Aleutians (Caricchi et al., 2014) and Santorini, Greece (Papageorgiou et al., 2019). Following this mechanism, the difference in depth between the sources responsible for the 2002–2005 uplift and 2009–2011 subsidence at Changbaishan does not imply the occurrence of two separated reservoirs but a single, vertically extended magma chamber affected by an early magma intrusion and pressurized degassing and by a later crystallization associated to passive gas release. Both the source depths retrieved for the 2002–2005 uplift and 2009–2011 subsidence are consistent with the 4–8 km deep basaltic trachyandesitic to trachytic magma chamber as discussed above for the 2018–2020. The depth of the point-source related to the subsidence is consistent with that determined for the shallower reservoir related to the 2018–2020 uplift (Supplementary Tables S1, S2).

Our results highlight that the Changbaishan magmatic system is still active and consists of reservoirs located at different depths whose present-day dynamics is associated with recharge processes, a feature consistent with the available petrological and geochemical data on past eruptions (Zhang et al., 2018; Yi et al., 2021) and geophysical data (Choi et al., 2013; Qiu et al., 2014; Hammond et al., 2020). A new upraise of magma may alter a part or the whole plumbing system and potentially trigger an eruption. The ALOS-2 data and related modeling results allow us to suppose the mechanisms of the present-day magma dynamics. Also, the successful exploitation of remote sensing data demonstrates their usefulness in case of under-monitored or partially monitored volcanoes. Because our data recorded a 2018–2020 deformation concentrated on the North Korea sector of the volcano, it may be that the in-situ monitoring system located on China territory does not detect the full pattern of signals related to the volcano dynamics. This prevents possible inferences and the correct interpretation of such signals, reducing the ability to properly evaluate the hazard and the efficiency of early warning systems.