Pyroclastic Density Current Hazards at the Baekdusan Volcano, Korea: Analyses of Several Scenarios from a Small-Case to the Worst-Case Colossal Eruption

The Baekdusan volcano was formed through three stages of activity: (a) a basalt shield (aging between 22.6 and 1.48 Ma), (b) a trachytic comendite stratocone (aging between 1.19 and 0.02 Ma), (c) a trachyte-comendite ignimbrite deposits (aging from 20 ka till date). Volcanic seismicity, ground deformation, and volcanic gas geochemistry yield new evidence for magmatic unrest of the volcano between 2002 and 2006. The monitoring data suggest that Mt. Baekdusan is a potentially active volcano and that its close monitoring is needed. One of the possible volcanic hazards from this volcano is the pyroclastic density currents. In order to evaluate the small-scale pyroclastic flow emplacement scenario of the 1903 AD eruption, Titan2d mass-flow model is used. The 1668–1702 AD and the Millennium eruption are characterized by 4–5 and ~7 VEI, respectively. The Millenni‐ um eruption can be considered as the last colossal super-eruption like Tambora, and so Baekdusan volcano could have even a global effect. The parameters used are as follows: volume (5–10 × 10 7 , 1 × 10 9 , 2 × 10 10 m 3 ), the vents of the 1903, 1668–1702, and Millenni‐ um eruptions; the pyroclastic flow runout calculated in the field are small (3,000 m : 1903 eruption), intermediate (5,000 m :1668–1702 eruptions), and large (7–80,000 m, Millennium eruption). The initial velocities (m/s) range from 50 (1903 eruption) to as high as 300 (Millennium eruption). The input parameters have constructed three scenarios (1903, 1668–1702, and Millennium eruptions) following the recent volcanic history of Baekdusan volcano. These trachyte and comendite lavas, and their pyroclastic equivalents during the cone-construction stage, and mostly comenditic pyroclastic rocks in the of unstable loose ig‐ nimbrite. The current may be a single-surge or more sustained where the collapse is retrogressive. Most large-volume ignimbrites derive from current types (B) and (C), which may involve periods of quasi-steady flow. Many may include significant components derived from currents of type (F) [44, 45]. For Baekdusan, we chose the following three types: Soufriere, Pelee, and Merapi types (modified from [62]).


Introduction
When active volcanoes are located in strategic and controversial territory with a history of political and military conflicts [1,2], their studies and the better understanding of how those volcanoes work are much more difficult, especially when the earth scientists have to deal with mitigation measures to protect the population living around [3][4][5][6]. This is surely the case of an active volcano known by a few until 15 years ago, situated between the China and North Korea borders, called by different names by the countries that surround it. The name of the volcanic field surrounding Tianchi(=Cheonji) caldera lake is 白頭山 in Chinese characters which means "White Head Mountain" These characters are Romanized as Baegdu-san, Baekdu-san, or Paektusan in Korea, Hakuto-san in Japan, and Changbaishan(長白山) in China [7, and references therein]. From now on, we will use the Korean name, Baekdusan, to avoid confusion. Baekdusan volcano has recently become of greater interest to volcanologists after the sign of unrest during the period 2002-2009, with several volcanic precursors such as volcanic earthquake, surface inflation, specific volcanic gas emission, higher temperatures of hot springs [8][9][10], which have prompted the Eastern Asian Volcanological Community to pay major attention to this volcano by establishing a Volcanological Observatory (the Chinese Tianchi Volcanological Observatory has collected data since 1999 [11]) and by setting up an improved monitoring network around the volcano [12,13]. Since 2012, the South Korean scientists have also taken steps in this direction, and the National Emergency Management Agency of Korea (Ministry of Public Safety and Security of Korea) has sponsored a vast research project on Baekdusan volcano (NEMA-BAEKDUSAN-2012-1-2), granted from the Volcanic Disaster Preparedness Centre (MPSS-NH-2015-81), and through the Natural Hazard Mitigation Research Group. Under the auspices of these groups, we analyzed the hazards from pyroclastic density currents and their impact on the population and infrastructures around Baekdusan volcano [14]. Further, we present arguments as to why this volcano has many characteristics as a potential site for the development of the conditions that can cause a super-eruption. Thus, we must also consider this extreme event if we do have to prepare for the next eruption, as it would have a global impact [15].

Regional tectonic of Baekdusan volcano and an outline of the plumbing system
Baekdusan volcano is an intraplate stratovolcano included in the larger Changbaishain area, where we can list a series of intraplate volcanoes (Wangtiane volcano, Bukpotaesan volcano, Zhenfengshan volcano, [7], Figure 1). Cenozoic intraplate basalts are widely distributed around the Changbaishan volcanic area. Beneath the Baekdusan volcano, a prominent low-velocity anomaly with a plume-like shape has been imaged in the upper mantle by P-wave tomography [16,17]. This suggests an upwelling of hot and wet material from the mantle transition zone (410-670 km [18]), where a stagnant subducted Pacific slab produces fluids and deep earthquakes through faults preserved in the stagnant slab. Hence, fluids and magmas may be supplied to the upper mantle under Baekdusan volcano, unlike beneath other volcanic and nonvolcanic areas in NE Asia, because large deep earthquakes occur frequently in the vicinity (a distance of about 200 km) of the Baekdusan volcano. At least some of these earthquakes are caused by the reactivation of faults which were produced within the plate, even after plate subduction. Further, the fluids could be released to the overlying mantle wedge from the preserved faults in the slab when large deep earthquakes occur. As a result, the Baekdusan volcano has to be considered much more active and prominent than other intraplate volcanoes in NE Asia [9,10], because more magmas are produced from the fluids released from the slab by the deep earthquakes [19]. Further, Kuritani et al. [20] use the trace element and Pb isotope data to demonstrate the intense hydration of the big mantle wedge beneath the Baekdusan volcano, through two different dehydration events of an ancient stagnant slab (~1.5 Gyr) and the most recent subducted Pacific slab (Figure 2A and B show the most recent setting). The wet mantle wedge is likely to result in more frequent rise of diapiric wet melts, as shown from the tomography [17, references therein], and develop a deep hot zone intruded by mafic sills repeatedly at the MOHO depth or/and scattered in the lower crust [21, and references therein]. Enrichment in fluid-mobile elements (e.g. Cs, Ba, K, Sr, U, Pb), a lack of enrichment in nonfluid-mobile elements (e.g. Nb, Ta, Ti, etc), together with specific isotopic signatures (e.g. 87 Sr/ 86 Sr, 206 Pb/ 204 Pb) in island-arc basalt or subduction-related magmatism are interpreted as a result of metasomatism of the magma source by slab-derived fluids or melts. Their detailed  investigation has helped to establish the current, rather empirical, model of slab dehydration and subduction-related melting processes. But, according to petrochemical characteristics and magma genesis of Holocene volcanic rocks at Mt. Baekdusan [22], trace elements from the Baekdusan volcano do not indicate significant contributions from a subduction slab of the Pacific plate. Volcanic activity in Mt. Baekdusan summit during Holocene resulted from intraplate magmatism [after 23,24], caused by upwelling of asthenospheric mantle under extensional tectonic conditions (Figure 3 [24]). Hot zone melts are H 2 O-rich as in Baekdusan volcano (e.g., Millennium eruption [25]). Consequently, they have low viscosity and density, and can readily detach from their source and ascend rapidly. Crystallization begins only when the ascending magma intersects its H 2 O-saturated liquidus at shallow depths. Decompression and degassing are the driving forces behind crystallization, which take place at shallow depths  [17], DXAL, Daxing-Anling (modified from [17]).
Updates in Volcanology -From Volcano Modelling to Volcano Geology on timescales of decades or less. Degassing and crystallization at shallow depths lead to large increases in viscosity and stalling of the magma to form volcano-feeding magma chambers [26,27]. The above model is well confirmed from the study of Wei [24] at the Baekdusan volcano. The heating of the deep crust was inefficient in early times, because much of the basalts erupted out directly, and much of the heat from magma was not used to heat the wall rocks. Annan and Sparks [27] state that there is an initial incubation period in which the basalt intrusions solidify. Generation of silicic melts initiates when the solidus temperature of either the basalt magma or the surrounding crust is reached. At an intrusion rate of 50 m per 10,000 years, incubation periods in the range 10 5 -10 6 years are estimated, consistent with geochronology and stratigraphy at Baekdusan volcano, with its evolution from mafic to silicic volcanism (the stratigraphy and geochronology will be treated in the next section). Once a zone of residual partial melt develops in the deep hot zone and starts to prevent the basalt from eruption, the temperature rises more efficiently, and the production of residual melt is accelerated. Large amounts of evolved residual silicic melt have been formed in the last 10,000 years and then intruded into the upper crust to form magma chambers large enough to result in the calderaforming explosive eruptions. The repose time between eruptions of size similar to the Millennium eruption was about thousands of years, one order greater than less-explosive parasitic eruptions of hundreds of years [7].

Summary of stratigraphy and geochronology of Baekdusan volcano
The volcanic activity at the Baekdusan volcano goes back as far as nearly 25 Ma B.P [7,28,29]. Wei et al. [7] describe the volcanic activity in three stages: early shield-forming; middle coneconstruction; and late ignimbrite-forming eruptions. The volcanic rocks were alkaline and tholeiitic basalts during the shield stage, trachyte and comendite lavas, and their pyroclastic equivalents during the cone-construction stage, and mostly comenditic pyroclastic rocks in the final stage (Table 1, Figure 4). However, Yang et al. [30] affirm on the basis of a large amount of age and stratigraphic data ( 40 Ar/ 39 Ar, U-series disequilibria) that Baekdusan volcano has been more active in the last 20 ka (~19, 16, 11 ka) than that was previously thought. From historical records, several authors [31,32] affirm the high frequency of eruptive events, and since the Millennium eruption, more than 31 eruptive events have been documented, most of which are the Plinian eruptions with volcanic ash that dispersed into the regions in the vicinity of the volcano. Many authors concord that the last explosive eruptions are in chronological order, 1668 AD ( Figure 5), 1702 AD (Figure 6), and 1903 AD (Figure 7) [32][33][34]. In Figure 4 is shown the Baekdusan volcano's geological map with all the terrains outcropping (see the legend of Figure 4 for more major details). The shield-forming stage lasted from the Miocene into the early Pleistocene. Alkali olivine and tholeiite basalt form the main part of a shield-like lava plateau that covers an area of 7,200 km 2 centered at Cheonji. In addition to the eruptions from central vents, there also were fissure eruptions, as evidenced by feeder dikes along the headwaters of the Tumen (Duman) and Heishihe Rivers. The youngest shield basalts are located near Cheonji cone and are overlain by the alkaline trachytes of the cone, except in the northwestern part of the shield. The shield basalts are ∼150 m thick near the Baishan Forestry Center in the headwaters of the Yalu River (Aprok River). The cone-construction stage lasted from ∼1.0 Ma to perhaps as young as ∼20 ka. The latter is about the age of the youngest trachyte lava on the caldera rim near Tianwenfeng (K-Ar, 19 ± 5 ka) and one of the youngest comendite lavas, Qixiangzhan, which was dated recently at 17 ± 1 ka. This stage was characterized by eruption of trachyte and comendite lavas and pyroclastic rocks, and dominated by trachyte lava effusion from the center of Cheonji volcano. The cone-construction stage has been subdivided into four stratigraphic formations for the main cone and the Laofangzixiaoshan  [7,28,29].
Updates in Volcanology -From Volcano Modelling to Volcano Geology  part of the shield ( Table 1). In general, outcrops of lava decrease image with elevation on the cone. In plane view, the eastern sector of the cone consists of younger trachytes, whereas older lavas are distributed more widely. However, all the rocks share similar geochemical compositions and appear to be genetically related. The young Laohudong basalts are sparsely distributed as monogenetic vents around the caldera. Its largest and most explosive Plinian eruption occurred around 946-947 AD and is known as the "Millennium eruption."  This eruption has a volcanic explosivity index (VEI) of 7, comparable to the 1815 eruption of Tambora in Indonesia. Horn and Schmincke [25] have suggested that the Cheonji eruption column reached ca 25 km in height, injecting about 9 ×10 10 kg of volatile halogens into the atmosphere, with an estimated volume of pyroclastic material on the order of ~100 km 3 . Volcanic ash from the eruption has been found 1,000 km away from the site of the eruption. A 5 cm thick deposit was observed in southern Hokkaido, Japan [38]. Accordingly, the Millennium eruption of the Baekdusan Cheonji volcano has been considered as one of the greatest volcanic eruptions of the past 2000 years.

An outline of the Millennium eruption
One of the most largest destructive eruptions from the Baekdusan volcano occurred on Earth in the last millennium, about 1,000 years ago, and it is named the Millennium eruption [25,29,34,35]. It involved more than ~100 km 3 of magma erupted [10], and the Plinian column was estimated to have had a height of 25 km. A widespread ash-fall deposit has been traced across the Sea of Japan to Hokkaido, and it is several centimeters thick, at more than 1,000 km from the source [35,37]. Unwelded pyroclastic flows associated with collapse events of the Plinian eruption column extended to more than 70 km from the crater rim (Figure 8). The Millennium eruption formed a large caldera with a diameter of about 5 km and an area of 20 km 2 , forming the Cheonji (Tianchi = Sky) lake. Currently, the lake has a maximum depth of 374 m and an area of ~ 9.8 km 2 , with a water surface elevation of 2,189 m. It has a volume of 2 km 3 . A composite stratigraphic section of the Millennium eruption is given by Machida et al. ( [35]; Figure 9). According to Machida et al. [35], the section is broadly composed of five levels (A, B, C, E, F; for major details, see Figure 9). Figure 10 shows some unwelded Millennium pyroclastic flow outcrops with relative thicknesses and sites. However, there is a controversy if the complete volcanological section of the Millennium eruption has to be interpreted as a single continuous eruption rather than one with two distinct pulses, with the initial phase anticipating the major phases of tens or hundreds of years [39]. The scientific literature reports a large range of ages for the Millennium eruption [40]. The results show a variety of periods ranging from approximately AD eighth to fourteenth centuries, which are the dates of the Balhae(= Bohai) and Goryeo dynasties [40, and references therein]. The last geochronological data (radiocarbon wiggle-match dating) on the Millennium eruption show a much more restricted range, respectively, of 938-939 AD [41], 946-947 AD [42], and 1024 AD [28].

Scenarios of the 1903 AD, 1668-1702 AD, and the Millennium eruptions
The scenarios envisaged here are based on the large amount of data that have been collected by the Korean research project, CATER2011-5210 [36], and by the new literature produced [25, 32-34, among others] in order to assign a selected VEI (Volcanic Explosivity Index [43]). The main parameters are summarized in Table 2 (VEI, maximum volume erupted, and column height) for each eruption. The volumes considered for just the pyroclastic density currents are obviously less than the total volume and are listed in Table 3. Figure 11 shows that between several volcanic-related disasters, the pyroclastic density currents have to be considered the most dangerous, given the number of fatalites in the last 2,000 years.  Table 2. Summary of the main parameters for the 1093 AD, 1702 AD, 1668 AD, and the Millennium eruptions. Volumes are taken from Newhall and Self [43] compared with the VEI for each eruption.
Branney and Kokelaar [44] illustrate in detail how different types of pyroclastic density currents originate (Figure 12). From Figure 12, the pyroclastic density currents most likely to occur at Beakdusan volcano are the following: A, B, C, depending on the size and the amount of interaction with external water. For laterally moving systems, the pyroclastic density currents can be identified as two end-member types: (A) Concentrated currents or pyroclastic flows, (B) Dilute currents or pyroclastic surges (Note that these are often a spectrum, with gradations in between). The end-member A has the following main characteristics: 1. It has solids in concentrations of tens of volume percentage, and so they have higher densities than surges.

2.
They have a free surface, above which solids' concentration decreases sharply.
3. Transport of material is by fluidization, and most flows are considered to be laminar.
Updates in Volcanology -From Volcano Modelling to Volcano Geology 4. Velocities vary, typically by 10s of meters per second up to several hundred meters per second, inferred from the heights of obstacles overcome by flows.
The end-member B has the following main characteristics: 1. It contains less than 0.1-1.0% by volume of solids, even near ground surface, and so they are relatively of low density.

2.
They are density-stratified, with the highest particle concentration near the ground surface.

3.
Material is transported primarily by turbulent suspension.

4.
Transport system is modeled as one that loses particles by sedimentation, which depletes the system of mass. Eventually, the system may become buoyant, in which case it becomes a plume.
In fluid dynamics, the equations that govern the concentrated current are well known [49,50], and the mathematics based on this volcanic phenomena is broadly used to model pyroclastic flows, along with other gravity currents: avalanches, seafloor turbidity currents, lahars, and lava flows. Two numerical codes apply these equations with small differences: Titan2D [51] and Volcflow [52,53], but so far they have not produced any modeling for diluted currents. The only numerical code that models the diluted currents is called PDAC [54]. An application is shown for Somma-Vesuvius volcano, where pyroclastic density currents are modeled on the data of the 1631 AD sub-Plinian eruption ( [55]; Figure 13). Kelfourn [53] tried to model twofluid model block-ash flows and cloud-ash surge at Merapi volcano, developing new equations, but there is not yet an agreement (Kelfourn, personal communication). The importance of modeling diluted currents together with dense currents for the hazard of an explosive active volcano is commonly witnessed from the volcanological history, as seen in the work of Mastrolorenzo [56] who demonstrate the causes of mortality in PDCs at Pompeii and surroundings, on the basis of a multidisciplinary volcanological and bioanthropological study. Field and laboratory studies of the eruption products and victims, merged with numerical simulations and experiments, indicate that heat was the main cause of death. Their results show that exposure to at least 250°C hot surges at a distance of 10 km from the vent was sufficient to cause instant death, even if people were sheltered within buildings. Despite the fact that impact force and exposure time to dusty gas declined toward PDCs' periphery up to the survival conditions, lethal temperatures were maintained up to the PDCs' extreme depositional limits. This study suggests that diluted currents are very dangerous, and if the (F) Deposit-derived pyroclastic density current caused by gravitational collapse and avalanching of unstable loose ignimbrite. The current may be a single-surge or more sustained where the collapse is retrogressive. Most large-volume ignimbrites derive from current types (B) and (C), which may involve periods of quasi-steady flow. Many may include significant components derived from currents of type (F) [44,45]. For Baekdusan, we chose the following three types: Soufriere, Pelee, and Merapi types (modified from [62]).
Updates in Volcanology -From Volcano Modelling to Volcano Geology runout of the pyroclastic flow must be chosen, it is always underestimated when considering diluted currents, especially for large eruption.

Relationship between the Baekdusan eruptions: 1903 AD, 1668-1702 AD, and Millennium eruptions, and the most devastating eruptions in the recent history of the earth
As indicated in Table 1, the 1903 AD, 1668-1702 AD, and the Millennium eruptions are characterized by the following VEI, respectively: 2-3, 4-5, and ~7. From the classification of Newhall and Self [43], the magnitude of each eruption with all the volcanological parameters can be deduced. If we compare the total volume ejected from Baekdusan volcano in the last Millennium activity with the most devastating eruptions in the recent history of volcanology, it appears that the two most recent eruptions shown, from Mount St. Helens and Mount Pinatubo, each caused extensive damage and received intense media interest. Ejected volumes, however, were only about 1/10,000 (St. Helens) and 1/500 (Pinatubo) of the volumes associated with either the 74 ka Toba or the 2 Ma Yellowstone super-eruptions. Even the 1883 Krakatau eruption, which accounted for over 35,000 deaths, was smaller than the Toba and Yellowstone eruptions by two orders of magnitude. The 1815 Tambora eruption, which had an ejected volume of less than 5% of those of Toba or Yellowstone and was distinctly less than "super," affected global climate and was responsible for over 50,000 deaths ( Figure 14A and B [57]).  Figure 14B). The Millennium eruption may have been slightly larger than the 1815 Tambora eruption, Indonesia. Millennium eruption can be considered a super-eruption alike other eruptions (e.g., Campanian Ignimbrite, Campi Flegrei, Napoli [58]; Figure 15). This assumption is confirmed from some South Korean scientists that during the 2nd International Workshop for Volcanic Disaster Preparedness at Chungbuk National University, Korea affirmed that the Millennium eruption could have been as high as 7 VEI. Further, Xu et al. [42,59] suggest that the Millennium eruption was a Toba-like "ash giant/sulfur dwarf" and had much smaller global climatic impacts. Sun et al. [60] confirm that the eruption most probably occurred around AD 940s, 7 years after the Eldgjá eruption on Iceland, and examining the eruption's potential for climate-forcing using the sulfate records Updates in Volcanology -From Volcano Modelling to Volcano Geology from the ice cores, they conclude that it was unlikely to have had a global impact. Zircons U-Th data of Millennium eruption indicate that magma resided only ~8 ka prior to eruption. Therefore, Baekdusan volcano can produce catastrophic, explosive eruptions in the foreseeable future [36]. Thus, Baekdusan is the most active volcano in China [14], and it is a high-risk volcano because more than 100,000 people live on or near the slopes, with the addition of many tourists.

An example of Titan2D simulation with the size of 1903 AD eruption
The size of the 1903 AD eruption has been evaluated, and consequently some Titan2D simulations are presented, taking into account where the vents of the 1903 AD, 1668 AD, and 1702 AD were located. Five vent locations are considered: three inside the caldera lake and two outside the caldera, as illustrated in Figure 16. The simulations start with an initial column (or pile) that collapses. The column is a vertical cylinder with a diameter 250 m and the height Hcol ( Table 4). In Figure 17 we report the thickness of the deposit when the flow is stopped. We consider 1 h (3,600s) after the initial collapse of the eruption column which can be formed by the Plinian eruption. The scenario performed shows, according to the vent locations, that the flow moves in diverse directions (NE, SE) with a thickness of 3 m, and that if the vent is located inside the caldera the flow fills the caldera with a thickness of 5 m. If the vent is on the border of the caldera, the flow will deposit thickly, mostly in the northern valley, the upper stream region of Erdaobaihe (Figure 17 [34,61]).   Table 4. Parameters used to simulate the pyroclastic flow of 1903 AD eruption needed by Titan2D code. Figure 17. Results of the Titan2D simulation for a small-scale scenario (VEI 3) similar to the 1903 AD eruption [34,61].
Updates in Volcanology -From Volcano Modelling to Volcano Geology

The impact of 1668-1702 AD and Millennium eruption around Baekdusan volcano
From the pyroclastic density currents scenario of the 1668-1702 AD eruptions with a runout around 5,000 m (Figure 18), two densely populated cities which are nearest to Baekdusan volcano would suffer significant damage: Erdaobaihe and Songjianghe. Possibly, the pyroclastic density currents could also affect two other cities located slightly further away, such as Fusong and Lushuihe. In the case of the Millennium eruption, considering only the direct impact, without considering the indirect effect which in the case of such eruption could be severe, and considering an underestimation of the runout, the following cities will be hit: Erdaobaihe, Songjianghe, Changbai, Huishan, Baishan, Fusong, Lushuihe, and Tumer ( Figure  19).

Concluding remarks on the worst-case scenario
The worst-case scenario is surely an eruption of a colossal size (Millennium eruption), but it has to be excluded from the scenario of the 1668-1702 AD eruption. In this case, the hazard and the risk must be evaluated carefully, and the population around Baekdusan volcano must be informed in detail. We think in this case, the population at risk is not aware of the possible hazard of the volcano. In all sorts of circumstances, the population has to be alerted. Every sign of unrest has to be monitored and communicated to everybody involved in the decisionmaking to the people. As far as we are concerned there is not a robust monitoring system around Baekdusan volcano, especially on North Korea side. An emergency plan must be created for Baekdusan volcano.

Future research prospectives
In the near future, there are several points in research direction that the Asian scientists could undertake, if they have not done yet: (A) Evolution of the Baekdusan volcano plumbing system in the least 20,000 years. (B) Developing an open source numerical code to model diluted density currents (surge) for the worst-case scenario, (C) Making an event tree study of Baekdusan volcano like that of Somma-Vesuvius volcano, Napoli, Italy.