ARTICLE

https://doi.org/10.1038/s43247-023-01162-0 OPEN

Eruption plumes extended more than 30 km in

altitude in both phases of the Millennium eruption
of Paektu (Changbaishan) volcano
Antonio Costa 1 ✉, Leonardo Mingari2, Victoria C. Smith 3, Giovanni Macedonio 4, Danielle McLean3,
Arnau Folch2, Jeonghyun Lee5 & Sung-Hyo Yun5

The Millennium Eruption of Paektu volcano, on the border of China and North Korea, gen-
erated tephra deposits that extend >1000 km from the vent, making it one of the largest
1234567890():,;

eruptions in historical times. Based on observed thicknesses and compositions of the

deposits, the widespread tephra dispersal is attributed to two eruption phases fuelled by
chemically distinct magmas that produced both pyroclastic flows and fallout deposits. We
used an ensemble-based method with a dual step inversion, in combination with the FALL3D
atmospheric tephra transport model, to constrain these two different phases. The volume of
the two distinct phases has been calculated. The results indicate that about 3-16 km3 (with a
best estimate of 7.2 km3) and 4-20 km3 (with a best estimate of 9.3 km3) of magma were
erupted during the comendite and trachyte phases of the eruption, respectively. Eruption
rates of up to 4 × 108 kg/s generated plumes that extended 30-40 km up into the strato-
sphere during each phase.

1 Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Bologna, Bologna, Italy. 2 Geociencias Barcelona (GEO3-BCN), CSIC, Barcelona, Spain. 3 Research

Laboratory for Archaeology and the History of Art, School of Archaeology, University of Oxford, Oxford OX1 3TG, UK. 4 Istituto Nazionale di Geofisica e
Vulcanologia, Osservatorio Vesuviano, Naples, Italy. 5 Volcano Specialized Research Center, Pusan National University, Busan, South Korea.
✉email: antonio.costa@ingv.it

COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:6 | https://doi.org/10.1038/s43247-023-01162-0 | www.nature.com/commsenv 1

ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01162-0

L
arge explosive volcanic eruptions cover wide areas in thick eruption (ME) of Paektu, an active caldera situated on the border
pyroclastic deposits and release volatiles into the atmosphere of the Democratic People’s Republic of Korea (DPRK) and Peo-
that can perturb climate for years1–3. Although the impacts ple’s Republic of China (Fig. 1). The tephra erupted during this
of large volcanic eruptions can be far reaching, they vary con- event is found over large areas, including over Greenland that is
siderably. The impact on climate is linked to the amount of sulfur more than 7000 km from vent8, suggesting that it was a large
injected into the stratosphere, but models show that the latitude magnitude eruption of ~M6.7, but the climatic impact appears
and timing of the injection have implications for the longevity of limited9,10. This is at odds with what is observed from other large
the sulfate in the stratosphere and thus also control their eruptions (e.g., the 1815 CE eruption of Tambora, Indonesia11) so
effects2,4. Despite a wealth of detailed models, there are only a few it is not clear whether the eruption parameters for this event have
large eruptions for which there is enough data on the aspects been accurately constrained. Here we combine tephra fallout
thought to control climate impacts, and even still there are con- thickness measurements with a new ensemble-based method as
siderable gaps in understanding and large uncertainties5. Quan- part of a dual step inversion, which involves the FALL3D
tifying the volume of magma erupted, volatiles released and atmospheric tephra transport and deposition model12,13, to
plume height of a variety of past eruptions is critical to provide constrain the magnitude of the eruption and the eruption para-
empirical evidence to further interrogate the relationship between meters such as mass eruption rate, duration, and plume height.
explosive eruptions and their climate effects, as well as testing the For the first time we model the different phases of the eruption
accuracy of dispersal models. which are known from the proximal and distal deposits. By
Sulfate layers preserved in the polar ice cores provide a detailed coupling these new volume constraints with published petrolo-
chronology of volcanic eruptions in the last few thousand years, gical data on the volatile content14 we can update the estimates of
many of which have been attributed to particular eruptions using the amount of the volatiles released into the stratosphere and
glass compositions that act as a chemical fingerprint6. This time ultimately help understand the limited climatic impact.
anchored history allows for the impact of various eruptions to be
assessed using paleoenvironmental records (e.g., tree rings and A large eruption from Paektu. The discovery of a fine and
pollen7) and written records for more recent events. It is widespread volcanic ash (tephra) layer in northern Japan was the
important that these impacts are linked to eruption parameters to first evidence to suggest there was a recent eruption from Paektu
further inform our understanding on what controls the impact on volcano15, also referred to as Baekdu (South Korean), Baegdusan
climate. A notable eruption for which we have good constraints (Japanese), or Changbaishan (Chinese). This visible tephra was
on both the date and climate response is the 946 CE Millennium widely identified in sequences in northern Japan, including

Fig. 1 Studied area. Studied area and locations of tephra observations used for the dispersal modelling, compiled from terrestrial exposures and marine cores,
for the entire deposit and for the points where the two magmatic phases can be distinguished (see Supplementary Table 1 for further details and references).
The map was created using the QGIS open software (https://qgis.org/). Paektu caldera, around Lake Tianchi, in the inset is from astronaut photograph
ISS006-E−43366 acquired April 4, 2003 (ISS Crew Earth Observations experiment and the Image Science & Analysis Group, Johnson Space Center).

2 COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:6 | https://doi.org/10.1038/s43247-023-01162-0 | www.nature.com/commsenv

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01162-0 ARTICLE

southern Hokkaido, northern Honshu, and the Sea of Japan15. suggesting there may have been another eruption phase or other
Based on the decrease in tephra thickness to the east and the activity in early 947 CE.
presence of K-feldspar, consistent with observations from pyr-
oclastic density current (PDC) deposits around the volcano16,
Paektu was identified as the source by Machida and Arai15. This Proximal deposits of the Millennium eruption. The Millennium
distal tephra layer was labelled B-Tm following Japanese eruption stratigraphy around Paektu has been the focus of
tephrostratigraphy conventions, with B referring to Baegdusan numerous studies over the last three decades21. The studies have
(Paektu) and Tm relating to the type locality for the deposit at been focussed on the deposits on the separate sides of the border,
Tomakomai, Hokkaido15. Since these initial discoveries, the with seemingly no research group working on those from both
B-Tm tephra layer has been identified in numerous marine, China and DPRK. Various pyroclastic units are identified around
lacustrine and archaeological sequences across northern Japan, the volcano but their attribution to a particular eruptive event had
northeast China, and coastal regions of Russia (see McLean et remained somewhat ambiguous. Both Horn and Schmincke33
al.17 and references cited in Supplementary Table 1). In addition, and Pan et al.21 indicated the ME deposits had both a comendite
glass shards with B-Tm geochemical compositions have been and trachytic phase. However, Horn and Schmincke33 noted that
identified as a non-visible layer (cryptotephra) in ice cores from the comenditic component was the most dominant and only ~5%
Greenland8, over 7000 km from its source. The glass composi- of the melt erupted was trachytic in composition, found in min-
tions of the tephra deposits across Japan verify that the occur- gled clasts within the deposit. Their ratio of comendite to trachyte
rences are indeed from the same eruption of Paektu volcano and is inconsistent with the similar thickness measurements of the
have correlated the B-Tm tephra to the Millennium eruption comenditic and trachytic fallout units recorded by Pan et al.21,
(ME) deposits18. and at odds with the distal records that reveal that up to 50% of
the glass shards are trachytic in composition17. Pan et al.21
interpreted that all the proximal units labelled Baguamiao, Mil-
Geological history of Paektu volcano. Volcanism at the intra- lennium and some others, attributed by other researchers34,35 to
plate Mt. Paektu (42.00°N, 128.05°E, 2750 m a.s.l.) is fuelled by later eruptions, were all phases of the ME. As such, Pan et al.21
the upwelling of mantle material. Initial basaltic fissure eruptions labelled the deposits associated with the 946 CE eruption the
commenced at the location ~28 My ago and the volcanic edifice Generalized Millennium Eruption, here referred to as ME, and
started to develop at ~4 Ma with activity divided into three stages. note that the sequence includes a light grey and dark grey
A basaltic shield volcano developed in Stage 1 and more explosive coloured comenditic fallout and PDC and a much darker
activity began in Stage 2 at around 1 Ma and built the trachytic- coloured trachytic fallout and PDC. Both comenditic and tra-
rhyolitic stratovolcano19. The Cheonji (Tianchi in Chinese) chytic fallout phases are decimetres in thickness in sections
caldera-forming Stage 319,20 generated several large eruptions in ~30 km east of the vent. The comenditic PDC is more than 10 m
the last 80 ka21. Deposits are found in distal records testifying to thick at sections <15 km east of the vent and it is partially welded
multiple explosive eruptions during this most recent stage22–24, in places21,33. At least 0.4 m of a trachytic PDC was deposited
which have spanned a range of magnitudes from small vulcanian after a fallout deposit from the same magma, but it has been
eruptions, such as the 1903 event25, to the Plinian Millennium eroded in most sections suggesting that it is a minimum thickness
eruption26. for the deposit21. The observed deposit thicknesses both proxi-
mally and distally suggest that the comendite and trachyte fallout
phases were similar in size, but the comenditic PDC was greater
Precise date of the Millennium eruption. Radiocarbon dates of in volume relative to the trachytic PDC.
tree remains found with the B-Tm ash initially indicated the A lahar deposit was observed by Pan et al.21 between the two
eruption occurred between 794–1192 CE15. The date of the phases of the ME at a locality ~15 km east of the volcano. They
tephra was shown to be consistent with that of the deposits close suggested that this could reflect a hiatus of a few wet months,
to the source, which was subsequently constrained further over which is consistent with historic accounts that record phenomena
the years using radiocarbon measurements of organic material commonly related to volcanic eruptions from 3rd November 946
recovered in the ME deposits around the vent. A precise wiggle- through to 7th February 947 CE21. Although the distal
matched radio-carbon age of 946 ± 3 CE was obtained from a tree sedimentary records seemingly indicate a single layer15, the
felled by the eruption that was recovered from the PDC on the sedimentation rates at the sites are not high enough to record
western flank of the volcano9. Further confirmation that the detectable sedimentation between deposits of the different phases
eruption occurred in 946 CE was established by Oppenheimer if they were indeed separated by months. In addition, current ice
et al.10 by obtaining radiocarbon measurements and identifying core sampling routines also mean that datasets are unlikely to
the 775 CE cosmogenic event27 across a similar tree stump show a clear separation in the different phases.
recovered from the same location within the deposits. The B-Tm Tephra fallout volumes were recently reassessed by Yang
tephra is also found in Greenland ice cores and a date of 939–940 et al.36 based on known proximal and distal thicknesses of the
CE was assigned to the eruption using the Greenland Ice Core deposit. They dealt with the entire deposit thickness (both
Chronology8. Subsequent research established the GICC05 ice comenditic and trachytic) and estimated that the fallout equated
core timescale was offset by −6 years28 so the ice core date is to 13.4 to 37.4 km3 Dense Rock Equivalent (DRE) of magma36.
consistent with that obtained using the tree stump. The PDC volumes associated with the eruption were estimated by
Written records from across the region document phenomena Horn and Schmincke33 to be between 5.3–7.6 km3 DRE (12.3 to
consistent with a large volcanic eruption in both 946 and 947 CE, 17.5 km3 bulk volume) based on assumptions of their extent in
with references such as the “roll of drums within the sky” in 946 China. Recently, Zhao et al.37 mapped the PDC deposits in China
CE (record from the Koryǒ dynasty, Korean Peninsula29,30), and estimated they accounted for ~ 3 km3 DRE (7 km3 bulk
“white ash falls as snow in the night” in Nara, Japan, on 3rd volume). These previous estimates were reconciled by Yang
November 946 CE10,30,31, and two records from 7th February 947 et al.36 to constrain the total extra-caldera PDC volume to
CE that noted “sound in the sky, like thunder”30,32. These historic ~4.1–5.7 km3 DRE (9.3–13.0 km3 bulk) and estimated up to
accounts and the tree and ice records indicate that the eruption ~2.1 km3 DRE (maximum of 4.9 km3 bulk tephra) of intra-
certainly started in 946 CE, with the further historical accounts caldera PDC. These estimates indicate that the bulk volume of the

COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:6 | https://doi.org/10.1038/s43247-023-01162-0 | www.nature.com/commsenv 3

ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01162-0

entire ME deposit is between 40.2 and 97.7 km3, which equates to We have 23 locations for which we know thicknesses for both
17.5 to 42.5 km3 DRE magma assuming a tephra deposit density the comenditic (first) and trachytic phases of the ME
of 1000 kg/m3 and a magma density of 2300 kg/m3 that is deposit17,21,40. These data points were used for the inversion
appropriate for rhyolite. Nonetheless, these represent two (details below), and the solutions related to each eruptive phase
different magma compositions and the ratio of comendite to then used to compare and validate the results with the total
trachytic magma was not quantified. thicknesses of the deposit for which more points are available
(n = 83; see Fig. 1 and Supplementary Table 1).
Impact. Proxy records of past climate such as the maximum
Latewood Density of tree rings, which is sensitive to summer Solving an inverse problem using an ensemble-based method.
temperatures7, from a range of Northern Hemisphere locations, To estimate Eruption Source Parameters (ESPs) from observed
including China, do not show noticeable changes following the deposit data we used a dual step inversion method. The first-step
Millennium eruption10. Furthermore, the non-sea salt sulfur inversion used a fast semi-analytical code, such as PARFIT-2.3.1
(nssS; i.e. that deposited from volcanic emissions) recorded in the (http://datasim.ov.ingv.it/download/parfit/download-parfit-2.3.1.
Greenland ice core is ~30 ppb at 946 CE, roughly double the html), to identify suitable wind conditions and ranges of ESPs for
background levels, with no perturbation in nssS in Antarctic the subsequent ensemble-based inversion. For this purpose, we
records28. These nssS are less than half that recorded in Green- considered more than 18,000 daily-averaged wind profiles near
land (80 ppb; 80–130 ppb in Antarctic ice cores) following the the eruptive vent location (Supplementary Table 1) over a period
1815 CE eruption of Tambora, Indonesia, but the latitudes are of 50 years (1961 to 2010), which were obtained from the ERA5
very different which would have affected sulfate deposition28. Reanalysis41.
This muted climatic response indicates either low S volatile We solved three separate inverse problems, one for the
concentrations or tropospheric plume heights. comenditic phase, one for the trachytic phase, and another for
the total deposit (two phases together; Supplementary Table 1).
New insights into the Millennium eruption. Here we use all the For each case we used the weighted least square method,
published thickness and compositional data of the deposits that considering both proportional and statistical weights (i.e. θ1;i /
are found over an area >700,000 km2 east of the volcano (see 1=T i 2 and θ2;i / 1=T i , being θ the weights and T deposit loads42),
Fig. 1 and Supplementary Table 1) with a dual-step inversion and allowed wide ranges of ESP and meteorological conditions. In
method to model the tephra dispersal from the two phases of the this step, we considered only the tephra load at distances >40 km
Millennium eruption. A new ensemble-based inversion modelling from the vent where the comendite and trachyte fraction have
strategy38 was used to gain insights into the eruption dynamics been estimated (see Supplementary Table 1; proximal outcrops
(e.g., mass eruption rate and column height) and estimate the were excluded as the model is oversimplified at these distances43).
volume of the comenditic and trachytic magma erupted during The results of this step generated the ranges of input needed for
the ME. We then use these total volume estimations with existing the second step (see “Methods”).
data on the volatile compositions of the melts14,33 to estimate the The results of the first step inversion are reported in Table 1.
likely total volatile release. The probabilistic nature of the For the comenditic phase, the total erupted mass (TEM),
ensemble-based method38 used for solving this inverse problem expressed as DRE, ranges from 16 to 20 km3 and column height
allowed us to properly assess the uncertainties. is about 22–26 km above vent level (a.v.l.), whereas, for the
trachytic phase, the volume ranges from 17 to 36 km3 DRE with a
Results and discussion column height from 26 to 35 km (a.v.l.). Using the dataset
Tracing the dispersal of the comenditic and trachytic phases. consisting of all measurement points available at a distance from
Although the comenditic and trachytic phases of the ME deposit the volcano greater than ~40 km, we found an erupted magma
are visually distinct in proximal locations within 40 km of the volume of about 15–19 km3 DRE and a column height ranging
volcano, at distal locations (more than a few hundred kilometres from 38 to 45 km. The wide ranges obtained for the TEM and
from the vent) the two phases cannot be distinguished visually15. column height, and the fact that the total volume is smaller than
However, since the melts are compositionally distinct the glass each individual phases, reflects the intrinsic uncertainties in the
chemistry can be used to establish which phase is present at a site. solution of the inverse problem with simplified analytical models
Many distal sites have deposits from both phases, and the per- such as PARFIT, which are not well constrained due to the
centage of which phase is present can be estimated based on the scarcity and sparsity of the data44,45. Furthermore, solving the
proportion of glass shards of each composition17. We used this column height and wind intensity with only thickness data is not
published glass chemistry to attribute a thickness to each phase at ideal46,47. For these reasons, these first-step ranges were used as a
distal sites where the tephra thickness is also reported. This starting point to find an inverse problem solution using the
assumes that the glass compositions analysed are representative, second ensemble-based method.
which may not be the case, but given that they are similarly easy In the second step, for each optimal meteorological configura-
to analyse there is no likely bias. The B-Tm ash in some distal tion (Table 1) we performed six ensemble runs using the
records does not form a visible layer, with the layers (termed FALL3D-8.2 dispersal model38,48 considering an ensemble size of
cryptotephra) being identified through density extraction 128 members. The best meteorological configurations (start date
techniques39. McLean et al.17,40 carried out detailed cryptotephra reported in Table 1) were used to drive the dispersal model for
analysis of the Lake Suigetsu sediment cores (central Japan, each ensemble run. This approach provides an ensemble of
~950 km from Mt Paektu) and noted that the B-Tm is only visible deposit mass loading field estimations (first guess). Subsequently,
in some core sections and has 1 mm in thickness. The tephra observations are assimilated to reconstruct the deposit mass
corresponds to shard concentrations of 25,000 shards/gram dry loading (analysis) using the GNC method described by Mingari
sediment in the 1 cm (core depth) sample40. Based on these et al.48. The deposit was reconstructed for the individual phases
measurements, we have assumed that 25,000 shards/gram of dry considering the six meteorological configurations and assimilat-
sediment equates to 1 mm in thickness and used this to assess ing the 23 observations (assimilation dataset) corresponding to
likely tephra thickness where only shard concentrations are each phase. The total tephra deposit is obtained as the sum of the
known (see Supplementary Table 1). comenditic and trachytic deposits. The total deposit resulting

4 COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:6 | https://doi.org/10.1038/s43247-023-01162-0 | www.nature.com/commsenv

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01162-0 ARTICLE

Table 1 Results of PARFIT tephra dispersal inversion.

Eruptive parameters Explored range Combined Combined Comendite Comendite Trachyte Trachyte
(θ1 ) (θ2 ) phase (θ1 ) phase (θ2 ) phase (θ1 ) phase (θ2 )
Optimal ERA5 reanalysis 30 30 29 8 23 21
meteorological 1961–2010 September October November October September April
configuration 1968 1992 1978 1975 1986 1994
TEM Calculated 4.5 3.6 3.9 4.9 4.1 8.7
(1013 kg) analytically
Volume Calculated 18.9 15 16.4 20.4 17.4 36.5
(DRE; km3) analytically
Column height 15–45 38 45 26 22 26 35
a.v.l. (H, km)
Suzuki coefficient70 1–9 6 2 4 2 6 3
Diffusion coefficient 1000–30,000 8000 16,000 20,000 28,000 19,000 20,000
(K, m2 s−1)
Density of aggregates 100–600 300 300 300 300 275 300
(kg m−3)
Deposit density Assumed 1000 1000 1000 1000 1000 1000
(kg m−3)
χ2 Calculated 0.72 0.37 0.31 0.10 0.13 0.05

Results of PARFIT tephra dispersal (first step) inversion for the ME eruption. θ1 and θ2 represent proportional and statistical weights in the PARFIT inversion, respectively.

from every combination (36 possible configurations) was first-guess emission source term for the ith ensemble member (si )
evaluated using different metrics computed using the full dataset is then rescaled according to si ! wi si . The resulting total source
of 83 measurements of total deposit thickness (independent term (in kg s−1 m−3), defined as the weighted sum of the
validation dataset). individual source terms:
The root-mean-square error (RMSE) was used to measure the
differences between the values observed and the analyses s ¼ ∑i wi si ð1Þ
interpolated at the sampling sites. For this purpose, we used the is represented in Fig. 4 in terms of the linear source emission
weighted version of RMSE and other statistical indicators, as strength δs (in kg s−1 m−1), i.e. δsdA being dA the area of the cell
explained in the Methods section. An overview of all used grid. The emission rate (i.e. δs vertically integrated) is also
statistical parameters is reported in Supplementary Table 2. represented in Fig. 4 (dash-dotted blue line) along with the
The optimal configuration, that minimises the Confidence cumulative erupted volume (solid black line). The total erupted
Ratio and Composite index (that is on average all the other volume for the fallout phases according to the inverted source
metrics), corresponds to the deposit reconstructed using the term was ~7.2 and 9.3 km3 DRE for the comendite and trachyte
meteorological datasets starting on 23 September 1986 for the phase, respectively (see Fig. 4 and Table 2), with a mean error of
comendite phase and 30 October 1992 for the trachyte phase (see about a factor 2.4, that is a volume between 3 to 17 km3 for the
Supplementary Table 2). Similarly optimal solutions for the are comendite phase and between 4 and 22 km3 for the trachyte
comendite phase were obtained with the meteorological condi- phase.
tions starting on 23 September 1986, and additional optimal
metrological conditions for the trachytic phase were those
Volume and eruption parameters. The results from this new
observed on 21 April 1994 (RMSE and MAE) and 29 November
approach indicate that fallout (Plinian and co-ignimbrite) from
1978 (Bias and SMAPE) (see Supplementary Figures). In our
each of the compositionally distinct fallout phases erupted
approach these represent the most compatible meteorological
roughly the same volume of magma, that is about 7–9 km3 DRE
conditions scanned across the ERA5 dataset. The reconstructed
(Table 2). This indicates that around 16 km3 of magma was
deposit fields for the comendite and trachyte phases and for the
dispersed by the two fallout phases, but, considering the uncer-
total deposit are shown in Fig. 2.
tainty, such a value could have been as low as 7 km3 and up to
Figure 3 shows measurement and analysis data from the
39 km3. For comparison, Horn and Schmincke33 calculated a
assimilation datasets (23 measurements for each eruptive phase)
volume of 23 ± 5 km3 DRE for the comendite phase, and recent
and for the validation dataset (83 measurements for the total
re-estimates by Yang et al.36 for the fallout for both phases were
deposit). Note that a RMSE close to 1, like in the presented case
estimated to be between 13.4 and 37.4 km3 DRE. The tephra
(Supplementary Table 2), means that deviations between
fallout volumes need to be added to the PDC estimates, which are
measurements

and

analysis are similar to the observation errors, 4.1–5.7 km3 DRE for those deposited outside the caldera and an
i.e.
yobs
i  yi
an

 ϵi . On the other hand, a positive bias additional ~2.1 km3 DRE for PDC emplaced in the caldera36 to
(0.39–0.40) for the optimal configuration means that the analysis estimate the total volume erupted. Our new fallout volumes
tends to underestimate the observations, but the bias is within the combined with those estimated for the PDCs indicate the total
observation error range (jBiasj<1). Finally, more than 68% of the volume of magma erupted during the ME was about 23 km3 DRE,
observations are within a factor 2.4 (that is one standard with at least 13 km3 and up to 47 km3 DRE erupted, corre-
deviation for a lognormal distribution). The other two optimal sponding to a M6.5 to M7.0 eruption. The average volume
solutions are reported in the Supplementary Figures. (23 km3) is similar to the amount of material that missing from
The source term was inverted using the procedure based on the edifice (22 km3) that would have collapsed into the void gener-
GNC method described by Mingari et al.48. This method provides ated by evacuating the magma. The amount of material displaced
a set of weight factors (wi ) for every ensemble member (see was estimated using 2.8 km as the mean radius of the caldera,
“Methods” section). According to the inversion procedure, the obtained averaging the maximum distances along the different

COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:6 | https://doi.org/10.1038/s43247-023-01162-0 | www.nature.com/commsenv 5

ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01162-0

directions, and ~900 m height based on the difference between

the maximum elevation and the bottom of the lake. The agree-
ment between these values suggests that the optimal estimations
are more realistic than the upper bounds.
Our results indicate that each of the two phases lasted less than
a day (~15 and 20 h; Fig. 4) and peak mass flow rates were
estimated to be about 4 × 108 kg/s, which are typical for Plinian
eruptions49. Similar values were estimated from observations of
the climatic phase of the 1991 Pinatubo eruption50,51. These high
mass flow rates would have fed umbrella clouds51 and sustained
columns that, for both phases, reached ~30–40 km (Fig. 4 and
Table 2), which are well above the tropopause height of ~12 km at
this latitude. This indicates that almost all of the eruptive mixture,
including volatiles released during the climatic phase of the
eruption, were injected into the stratosphere.

Volatile release and impact. The petrological method52,53 was

previously used with the comendite volumes33 (23 ± 5 km3 DRE)
with volatile contents of the melt nclusions (MI) and matrix
glasses to estimate the amount of S, Cl and F released. Mea-
surements of the volatiles in MI and matrix glasses from the
comenditic magma have been made by Horn and Schmincke33,
Iacovino et al.14 and very recently re-analysed by Scaillet and
Oppenheimer52.The F and Cl contents of MI and matrix glasses
cover a similar range14,33, suggesting the melts were probably not
saturated in F or Cl, and loss of these volatile phases could be
negligible. Iacovino et al.14 modified the petrological approach
used to also estimate the syn-eruptive degassing to quantify the
total pre-eruptive volatile loss associated with the comenditic
magma. This approached assumed the comendite erupted was
linked to the trachyte erupted via fractional crystallisation, which
is consistent with regression modelling of the major elements.
They compared the observed volatile contents relative to
incompatible elements (U in this case) to values estimated from
fractional crystallisation models to estimate the volatiles lost from
the melt during crystallisation14. Their data indicates that the pre-
eruptive volatile release from the comenditic magma could have
been substantial.
If we consider comendite magma volumes calculated with the
new method presented here, which range from 3 to 17 km3 DRE,
and the mean volatile contents measured by Iacovino et al.14 in
the MI and matrix glass, we can rescale the volatiles that were
removed from the melt prior to eruption to be between 5 and 30
Tg S, 6–32 Tg F, and 2–15 Tg Cl. While our scaled estimates for
the syn-eruptive S released into the stratosphere is 0.5–2.6 Tg
(1.0–5.2 Tg SO2). Nonetheless, these estimates do not consider
the volatile release associated with the trachytic melt (4–22 km3
DRE) as there are no published volatile data for trachytic matrix
glass or MI in crystals exclusively associated with the trachytic
magma. The new volumes suggest the pre-eruptive volatile release
would have been smaller than previous reported by Iacovino
et al.14. The volatile elements released prior to the eruption are
likely to have been sequestered by the aquifers and hydrothermal
fluids or passively degassed, and it is expected that most of these
volatiles did not make it into the stratosphere. Nonetheless, it
Fig. 2 Reconstructed tephra deposits. Reconstructed deposit thickness appears that some of the pre-eruptive S loss from the 1991
contours for the comendite (a), trachyte (b) eruptive phases, and for the total Pinatubo magma made its way into the stratosphere (which is
deposit (c) of the Millennium eruption (ME). The maps were created using higher at the more tropical latitude; ~18 km) as the satellite
the open software “Cartopy: a cartographic python library with a matplotlib measured S content in the stratosphere was greater than the syn-
interface” (https://scitools.org.uk/cartopy/docs/latest/index.html). eruptive load estimated by the petrological method53.
The lack of recorded climatic impacts for the eruption7,10
suggests a low syn-eruptive volatile release, which is consistent
with the above scaled syn-eruptive yields despite not including
any release from the trachytic melt. The S yield is similar to the
~2 Tg estimated from the ice core non-sea salt sulfate record54,

6 COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:6 | https://doi.org/10.1038/s43247-023-01162-0 | www.nature.com/commsenv

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01162-0 ARTICLE

Fig. 3 Comparison models and observations. Quantitative comparison between observed and modelled (analysis) deposit thickness of the comendite (a)
and trachyte (b) eruptive phases of the Millennium eruption (ME) using the assimilation dataset (23 measurements for each eruptive phase). In addition,
the comparison between the total deposit thickness and the validation dataset (83 measurements) is also shown (c). Note that zero-valued observations
for the comendite (P82 and P83, Supplementary Table 1) and trachyte (P34, P68, P72, P74, Supplementary Table 1) phases are not shown in these log-log
plots (a) and (b), respectively. The plot corresponding to the total deposit (c) includes the full dataset of 83 data points. Errors are estimated in accord to
Mingari et al.38 assuming a minimum value of 30%. The plots were created using the open software Matplotlib (https://matplotlib.org/stable/).

COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:6 | https://doi.org/10.1038/s43247-023-01162-0 | www.nature.com/commsenv 7

ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01162-0

Fig. 4 Time evolution of the eruption columns. Time evolution of the emission source profiles derived for the comenditic first phase (a) and trachytic
phase (b) of the Millennium eruption (ME). The plots were created using the open software Matplotlib (https://matplotlib.org/stable/).

which represents both the comenditic and trachytic phases. The Despite the lack of any significant climate impact, there would
modelled climate impacts of high latitude eruptions with S loads have been quite a substantial regional impact linked to the
of <5 Tg are shown to be negligible, especially those in winter as deposition of the thick deposits, lahars, and gas release. It is likely
insolation affects the removal of swimschool@abingdon.org.uk- that the eruption led to acid rain, which would have affected
sulfate from the stratosphere55. Nonetheless, recent simulations crops, livestock, and high F potentially poisoned local freshwater
have shown that stratospheric sulfur injections in the extratropics supplies in the downwind regions.
spend longer in the stratosphere than those from low latitudes,
and suggest that cooling is more pronounced, but confined to the
hemisphere in which it is injected4. Stratospheric halogens also Conclusions
exacerbate the radiation forcing of the sulfate56, but measured We applied a dual step strategy for solving an inverse problem
values suggest the syn-eruption halogen release during the ME aimed at reproducing the observed tephra deposit thicknesses
was low. In addition, the halogens are probably scavenged in the generated by the ME of Paektu volcano. The first step uses a fast
lower plume with only a few percent making it into the efficient code for tephra sedimentation based on an analytical
stratosphere57. solution (PARFIT). The second step considers a recently

8 COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:6 | https://doi.org/10.1038/s43247-023-01162-0 | www.nature.com/commsenv

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01162-0 ARTICLE

Table 2 Eruption parameter model output. based on ensemble approaches in order to produce improved
simulations of volcanic clouds and tephra deposits consistent
with real observations.
Eruption parameters Comendite phase Trachyte phase
Total erupted volume (km3DRE) 7.2 (3–17) 9.3 (4–22) Inversion method. We used a novel strategy to model the tephra
Plume height (km) 30–40 30–40 dispersal of the Millennium Eruption that improves, from a
Peak mass flow rate (108 kg/s) ~4 ~4
computational point of view, the approach of Costa et al.59, which
Duration (h) ~15 ~20
was based on a set of several hundreds of simulations of the
FALL3D ash dispersal model12 and 3D time-dependent meteor-
proposed ensemble-based inversion method with an High Per- ological fields across the region of interest. This new approach
formance Computing (HPC) numerical tephra dispersal model involves two separate inversion steps. The first step uses a simpler
(FALL3D) to reconstruct the eruption dynamics and the volume analytical solution60,61 (http://datasim.ov.ingv.it/download/parfit/
of magma erupted during the comendite and trachyte phases. download-parfit-2.3.1.html) to identify compatible meteor-
Such estimates of past eruption volumes and parameters, and ological conditions that describe the general observed tephra
their uncertainties, are key to understanding the impacts on cli- dispersal and estimate the most likely ranges of the Eruption
mate. The new inversion method allows us, for the first time, to Source Parameters (ESP). These first inversion estimates are then
estimate the evolution of the mass eruption rate with time during used in the subsequent step built on the ensemble-based inversion
the two main phases, and not only the time-averaged values, as algorithm recently proposed by Mingari et al.38 to further con-
commonly done. The total erupted mass of the two Plinian strain dispersal and the ESPs. This dual inversion step method is
phases are similar, with best estimates of ~7 km3 DRE for the computationally more efficient than that used by Costa et al.59
comendite phase and ~9 km3 DRE for the trachyte phase. Con- and thus, allows a very wide range of meteorological conditions
sidering previous estimations for PDC deposits and intra-caldera and ESP to be explored and allows for more extensive analysis.
filling, we can estimate a total magma volume of about ~23 km3
DRE for the eruption, which corresponds to a M6.7, consistent First step inversion. The first inversion is based on a semi-
with other recent estimations. The determination of the con- analytical solution of the mass conservation equation, which is
tribution of the two different magmas could be used to better not computationally intensive and allows tens of thousands of
estimate the amount of released sulfur once volatile data from the simulations to be explored on a simple desktop computer. We
trachyte phase (~9 km3 DRE) is available. The estimations for the used the most recent version of the PARFIT software package
comendite phase made with the proposed method are about one (PARFIT-2.3.1 code: http://datasim.ov.ingv.it/download/parfit/
third smaller than the reference values used in the literature. A download-parfit-2.3.1.html) to solve an inverse volcanological
simple scaling suggests between 0.5 to 2.8 Tg syn-eruptive sulfur problem, which consists of finding mean wind profile, column
release, which is in agreement with S loads estimated from ice- height, Total Erupted Mass (TEM), and Total Grain Size Dis-
core records, and consistent with negligible effects recorded in tribution (TGSD), using thickness and grain-size measurements
climate proxy records. This research highlights that in eruptions from available eruption deposit outcrops. PARFIT minimises the
that tap compositionally distinct melts, with different volatile difference between the simulated and the real deposit loadings.
histories, it is important to understand the temporal relationship The fitting is performed using a least-squares method which
between changes in melt chemistry and eruption dynamics so compares measured and calculated deposits thickness and grain-
that the likely impact of past eruptions can be accurately assessed. sizes. The function to be minimised42 is:
1 N  2
χ2 ¼ ∑ θi Y obs;i  Y mod;i ð2Þ
Methods N  p i¼1
Computational methodology
Tephra dispersal model. Tephra dispersal simulations were per- where θi are weighting factors, N is the number of measurements,
formed using the latest version (8.2) of the FALL3D-8.2 p is the number of free parameters, Y obs;i denote the observed
model12,38,48. FALL3D-8.2 is an Eulerian model for atmospheric ground loads (kg/m2) and Y mod;i the values predicted by the
transport and deposition based on the solution of the so-called model. The choice of the weighting factors, θi , depends upon the
the Advection-Diffusion-Sedimentation (ADS) equations12. This distribution of the errors42. To accurately constrain the input
code has been significantly refactored over the last years, resulting parameters, it is necessary to know thickness and grain-size
in major improvements to the scalability and performance12. In values for a large number, N, of locations, such as N ≫ p. Without
addition, the model physics and numerics have been revised in a large statistical set of observational data there is a risk of finding
this version12,58 to allow aerosol and radionuclide dispersal parameters that are incorrect, especially when considering limited
modelling. The updated model includes new coordinate mapping portions of the deposit61 that represent the tail of the distribution
options, a more efficient and less diffusive solving algorithm, and relative to a particle class44,46. In order to reduce the degrees of
better memory management and parallelisation strategy based on freedom, avoiding these problems46,47, the TGSD was estimated
a full 3-D domain decomposition, suitable for HPC applications. following Costa et al.62.
A set of case studies that use the model, including the simulation Tephra dispersal simulations are performed using the HAZ-
of the dispersal of long-range fine volcanic ash and volcanic SO2, MAP model that is embedded in PARFIT. The HAZMAP
and the deposition and dispersal of tephra fallout and radio- program simulates the tephra deposit generated by the sedimen-
nuclides, were recently presented by Prata et al.58. tation of volcanic particles from a sustained column60,61. The user
In addition, the capabilities of FALL3D have been extended by defines the ranges and resolution steps of the eruption
making it possible to perform an ensemble of model runs. parameters. This defines a multidimensional grid with combina-
Ensemble modelling allows model uncertainties, associated with tions of all parameters. PARFIT performs a search on the grid for
poorly constrained input parameters or errors in the underlying the minimum χ2 between the simulated deposit and the field data.
meteorological driver, to be characterised and quantified. Real wind profiles are used in the PARFIT-2.3.1 model that can
Furthermore, it is possible to explore different ensemble-based be assessed from an ensemble of daily winds over the investigated
data assimilation techniques or inversion modelling methods area. In this case, wind data were obtained from ERA5 Reanalysis

COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:6 | https://doi.org/10.1038/s43247-023-01162-0 | www.nature.com/commsenv 9

ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01162-0

Table 3 Eruption parameter ranges for PARFIT. Table 4 TGSD estimated according to Costa et al.62, and
accounting for ash aggregation following Cornell et al.68, for
Eruption parameter Range the PARFIT first step inversion.
Total erupted mass (TEM) Calculated
Plume height 15–45 km (Φ) Diameter (mm) Density (kg/m3) Shape Percentage
Suzuki coefficient70 1–9 (−5) 32.0 450.0 0.83 1.8
Diffusion coefficient 1000–30,000 m2/s (−4) 16.0 450.0 0.83 1.7
TGSD Following Costa et al.62 (−3) 8.0 450.0 0.83 2.2
Aggregate density 100–1000 kg/m3 (−2) 4.0 600.0 0.83 2.6
(−1) 2.0 800.0 0.83 2.6
(0) 1.0 900.0 0.83 2.5
dataset41 and converted to the required format using specific (1) 0.5 1100.0 0.83 3.1
routines within the software package. Due to the enhanced (2) 0.25 1300.0 0.83 6.4
computational efficiency, tephra dispersals driven by tens of (3) 0.125 1500.0 0.83 13.9
thousands of daily-averaged wind profiles were investigated to (4) 0.0625 1700.0 0.83 10.7
identify the most compatible meteorological conditions within (5) 0.03125 1900.0 0.83 5.3
the range explored (50 years, from 1961 to 2010). (2.5) 0.2 300.0 1.00 47.2
For this first-step inversion we considered thickness data
associated with the comenditic and trachytic phases, reported in
Supplementary Table 1, and searched for the best fit eruption function:
   T    T  
parameters for each distinct phase (similarly to Marti et al.63) in J y / y  y P1 y  y þ yO  y R1 yO  y ð4Þ
the ranges of values shown in Table 3. Since there is extremely
limited grain-size data, we estimated the TGSD using the Costa where the vectors y and y are defined as y ¼ Hx and y ¼ Hx,
et al.58 parameterisation that is based on magma viscosity and being H a linear observation operator which translates the model
column height. For this purpose, we considered a magma state into the observation space; P is the ensemble-based error
viscosity of about 107 Pa s for a rhyolitic composition14,17,19,64 covariance matrix computed in the observation space and we
and a plume height of 30 km, a typical value for Plinian assume a Gaussian distribution with error covariance matrix R for
eruptions49. We used the parameterisation proposed by Bona- the measurements. In this work, the observation operator
donna and Phillips65 to characterise densities of particles from interpolates the modelled mass loading into the sampling site
rhyolitic magma and used particle shapes that have been and computes the deposit thickness assuming a typical bulk
measured for other explosive eruptions66. Ash aggregation67 density for consolidated deposits44 of ρb ¼ 1000 kg/m3. This
was considered using the parameterisation of Cornell et al.68. optimisation problem leads to a problem of non-negative
The resulting TGSD used at this step is reported in Table 4. quadratic programming that can be solved using an iterative
approach. See Mingari et al.38 for further details.
Second-step ensemble-based inversion. As mentioned above, the The root-mean-square error (RMSE) was used to measure the
parameter ranges explored in this second step are derived from differences between the values observed and the analyses
the results obtained in the first step. The second step inversion interpolated at the sampling site:
reconstructed the spatial distribution of tephra deposit thickness rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1 n 2  obs 2
(i.e. mass loading) associated with the ME using the ensemble- RMSE ¼ ∑i¼1 pi yi  yan i
ð5Þ
based assimilation method proposed by Mingari et al.38. This n
algorithm provides a map of a tephra-fall deposit given a dataset being yobs
i the ith observation and yan i the corresponding analysis
of observations of a volcanic deposit at scattered points. To this for each individual phase. Notice that if the non-weighted
purpose, an ensemble of model realisations is additionally (pi ¼ 1) version of (1) is used, only very proximal data (i.e. largest
required. Given a background ensemble of model states char- deposit thickness measurements) makes a significant contribution
acterised by the state vectors xbi (i = 1,…, m), representing m to RMSE. Instead, in order to treat the deviations more evenly,
model realisations (ensemble members), along with a vector of the observation errors (ϵi ) are used to define the weights
observations y~o , the state estimate (analysis) is expressed as a according to pi ¼ 1=ϵi .
linear combinations of the background ensemble: Similarly, we computed the bias using the following definition:
1 n  
Bias ¼ ∑ pi yobs i  yi
an
ð6Þ
xa ¼ w1 xb1 þ ¼ þ wm xbm ð3Þ n i¼1
and the mean of absolute value of errors (MAE):
1 n


where the coefficients wi ≥ 0, i = 1,…, m are weight factors to be MAE ¼ ∑ pi
yobs i  yi
an

ð7Þ
determined which depend on the observations. In our case, the n i¼1
model state vector x is assumed to be the two-dimensional The mean absolute percentage error (MAPE) is also useful:
deposit mass loading on a regular grid. In the limit case in which


100 n
yobs i  yi
an

no observations are available, the ensemble members are uni- MAPE ¼ ∑
formly weighted according to the constant factors: wi ¼ m1 and, n i¼1 yobs
i
consequently, the state vector xa is simply the ensemble mean. where yobs is a vector of 83 observations and yan the vector of the
In order to determine the weight factors wi , a cost function corresponding analysis interpolated to the sampling sites.
must be minimised in the subspace spanned by the ensemble However, such definition of the error is very sensitive to the
members with the non-negative constraints wi ≥ 0. Let the vector observations having very small values (close to zero) and for this
model state ~x be a linear combination of the corresponding reason is preferable using the symmetric mean absolute
ensemble member vectors expressed as in Eq. (3), i.e. x ¼ ∑i wi xbi . percentage error (SMAPE), which is another accuracy measure
We look for the vector state minimising the following cost based on relative errors used to evaluate the performance of the

10 COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:6 | https://doi.org/10.1038/s43247-023-01162-0 | www.nature.com/commsenv

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01162-0 ARTICLE

Table 5 FALL3D model configuration. considered. Table 5 provides further details about the model
configuration.
A background ensemble with 128 members was generated by
Model parameter Value
perturbing Eruption Source Parameters (ESP), the horizontal
Simulation time 143 h wind components, and the aggregate density using either uniform
Ensemble size 128 or truncated normal distributions13,48. A Latin Hypercube
Grid resolution 0.23° × 0.16°
Sampling (LHS)13,69 was used to efficiently sample the parameter
Domain size 150 × 150 × 60
Emission source Suzuki distribution70
space. The ensemble configuration requires defining a central
shape value (i.e. results from the first inversion step), a perturbation
Species 10 ash bins + 1 aggregate bin range and a probability density function (PDF) for each
TGSD Estimated from magma composition and parameter to be perturbed. The ensemble is built by randomly
column height (Costa et al.62) sampling the corresponding PDF. Table 6 provides a summary of
the configuration used to define the background ensemble along
with the list of the perturbed parameters and model inputs.
Table 6 Ensemble configuration. The emission start time was uniformly sampled between t = 0
and t = 48 h, where t is the time after the simulation starts in
Model input Central value Perturbation PDF
hours, as detailed in Table 6. The emission source duration, 4t,
range was sampled in the range from 4 to 10 h, which means that
emission is only possible for 0 ≤ t ≤ 58 h. Since we aim to simulate
Column height 30 km avl 25% Uniform
Source start t = 24 ha 24 h Uniform
both phases of the ME (i.e. November 946 CE and February 947
Source duration 7h 3h Uniform CE events), we invert independently the two phases and consider
Suzuki 4 50% Uniform the formation of the total tephra fallout deposit as the sum of the
coefficient70 A first and second eruptive event for validation purposes.
MER Degruyterb – –
U Wind ERA5c 25% Gaussian Source term inversion. FALL3D solves an almost linear problem
V Wind ERA5c 25% Gaussian with weak nonlinearity effects (e.g., due to gravity current, wet
TGSD Costad – – deposition, or aggregation). Consequently, a re-scaling of the
Aggregate density 400 kg m−3 300 kg m−3 Uniform emission source term according to si ! wi si leads to a deposit
b
aTime in hours after the simulation start. mass loading re-scaled correspondingly: xbi ! wi xi , where si is
the emission source term in kg m−3 s−1 for the ith ensemble
bParameterisation from Degruyter and Bonadonna71.
cWind field obtained from the single level and model level ERA5 datasets.
dParameterisation from Costa et al.62, considering 10 ash bins from −3 Phi to 5 Phi and ash member. This is the mass injected into the computational domain
aggregates (200 micron diameter).
per unit time and depends on time and height above the vent. The
analysed tephra deposit (mass loading) on the ground is descri-
deposit reconstruction: bed by Eq. (3) that gives:


xa ¼ ∑i xbi ð10Þ
2 n
yobsi
y
an



i

SMAPE ¼ 100 ´ ∑
obs

an
ð8Þ i.e., the analysed mass loading is a result from the contribution of
n i¼1 yi þ yi
multiple members of an ensemble defined by the first-guess
emission source term si rescaled by wi .
We have also calculated the confidence ratio (CR), defined as
the value discriminating the symmetric band on the log-log space
of observations vs. analyses, containing 68% of the points
represented in Fig. 3 (this percentage will correspond to a sigma Data availability
The data used in this study are all published data and all the sources of information are
if we assume a lognormal distribution). Finally, we considered a reported in Supplementary Table 1 and cited in the Supplementary Reference section.
composite index (CI), defined as: Data sharing not applicable to this article as no original datasets were generated or
analysed during the current study.
CI ¼ RMSE þ Bias þ MAE þ CR þ SMAPE=100 ð9Þ

The values of the evaluation metrics described above are Code availability
reported in Supplementary Table 2. The Fall3D configuration files and scripts used for the simulations are available at https://
doi.org/10.5281/zenodo.10159965.
Ensemble modelling. The FALL3D computational domain was
configured using a horizontal resolution of 0.23° (longitude) and Received: 5 June 2023; Accepted: 5 December 2023;
0.16° (latitude) and a domain size of nx ´ ny ´ nz = 150 × 150 × 60
grid cells. Numerical simulations used the meteorological fields
estimated using the PARFIT-2.3.1 code. Among the different
meteorological configurations selected in the first step procedure,
the conditions corresponding to 25–29 January 1962, for the References
comendite phase, and 29 February–3 March 1976 for the trachyte 1. McCormick, M. P., Thomason, L. W. & Trepte, C. R. Atmospheric effects of
phase, showed the best agreement with observations during the the Mt-Pinatubo eruption. Nature 373, 399–404 (1995).
2. Robock, A. Volcanic eruptions and climate. Rev. Geophys. 38, 191–219 (2000).
second step (i.e. when using the ensemble-based FALL3D; see
3. Timmreck, C. Modeling the climatic effects of large explosive volcanic
Supplementary Table 2). The Total Grain Size Distribution eruptions. WIREs Clim. Change 3, 545–564 (2012).
(TGSD) was computed from the column height using the para- 4. Toohey, M. et al. Disproportionately strong climate forcing from extratropical
meterisation proposed by Costa et al.62 considering 10 ash bins explosive volcanic eruptions. Nat. Geosci. 12, 100–107 (2019).
with particle diameters ranging between −2Φ and 8Φ. Moreover, 5. Marshall, L. R. et al. Volcanic effects on climate: recent advances and future
avenues. Bull. Volcanol. 84, 54 (2022).
an additional bin aggregate of diameter dagg = 200 μm was

COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:6 | https://doi.org/10.1038/s43247-023-01162-0 | www.nature.com/commsenv 11

ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01162-0

6. Smith, V. C. et al. The magnitude and impact of the 431 CE Tierra Blanca 35. Yang, Q. F. The Distribution of Pyroclastic Flow on Changbaishan Tianchi
Joven eruption of Ilopango, El Salvador. Proc. Natl Acad. Sci. USA 117, Volcano. Project Report of Seismological Bureau of Jilin Province (in Chinese)
26061–26068 (2020). (2005).
7. Stoffel, M. et al. Estimates of volcanic-induced cooling in the Northern 36. Yang, Q. et al. The Millennium eruption of Changbaishan Tianchi Volcano is
Hemisphere over the past 1,500 years. Nat. Geosci. 8, 784–788 (2015). VEI 6, not 7. Bull. Volcanol. 83, 74 (2021).
8. Sun, C. et al. Ash from Changbaishan millennium eruption recorded in 37. Zhao, B., Xu, J., Yu, H. & Chen, Z. Pyroclastic density current facies in the
Greenland ice: implications for determining the eruption’s timing and impact. millennium eruption of Tianchi volcano northeast china: insights from
Geophys. Res. Lett. 41, 694–701 (2014). topography, stratigraphy, granulometry and petrography. Front. Earth Sci. 8,
9. Xu, J. D. et al. Climatic impact of the Millennium eruption of Changbaishan 72–91 (2020).
volcano in China: new insights from high-precision radiocarbon wigglematch 38. Mingari, L., Costa, A., Macedonio, G., & Folch, A. Reconstructing tephra fall
dating. Geophys. Res. Lett. 40, 54–59 (2013). deposits via ensemble-based data assimilation techniques. Geosci. Model Dev.
10. Oppenheimer, C. et al. Multi-proxy dating the Millennium eruption of 16, 3459–3478 (2023).
Changbaishan to late 946 CE. Quat. Sci. Rev. 158, 164–171 (2017). 39. Blockley, S. P. E. et al. A new and less destructive laboratory procedure for the
11. Oppenheimer, C. Climatic, environmental and human consequences of the physical separation of distal glass tephra shards from sediments. Quat. Sci.
largest known historic eruption: Tambora volcano (Indonesia) 1815. Prog. Rev. 24, 1952–1960 (2005).
Phys. Geogr. Earth Environ. 27, 230–259 (2003). 40. McLean, D. et al. Integrating the Holocene tephrostratigraphy for East Asia
12. Folch, A. et al. FALL3D-8.0: a computational model for atmospheric transport using a high‐resolution cryptotephra study from Lake Suigetsu (SG14 core),
and deposition of particles, aerosols and radionuclides—part 1: model physics Central Japan. Quat. Sci. Rev. 183, 36–58 (2018).
and numerics. Geosci. Model Dev. 13, 1431–1458 (2020). 41. Hersbach, H. et al. ERA5 hourly data on single levels from 1940 to present.
13. Folch, A., Mingari, L. & Prata, A. T. Ensemble-based forecast of volcanic Copernicus Climate Change Service (C3S) Climate Data Store (CDS). https://
clouds using FALL3D-8.1. Front. Earth Sci. 9, 74 (1841). doi.org/10.24381/cds.adbb2d47 (2018).
14. Iacovino, K. et al. Quantifying gas emissions from the “Millennium Eruption” 42. Costa, A. et al. Tephra fallout hazard assessment at Campi Flegrei caldera
of Paektu volcano, Democratic People’s Republic of Korea/China. Sci. Adv. 2, (Italy). Bull. Volcanol. 71, 259–273 (2009).
e1600913 (2016). 43. Armienti, P., Macedonio, G. & Pareschi, M. A numerical model for the
15. Machida, H. & Arai, F. Extensive ash falls in and around the Sea of Japan from simulation of tephra transport and deposition: applications to May 18, 1980
large late Quaternary eruptions. J. Volcanol. Geotherm. Res. 18, 151–164 Mt.St. Helens eruption. J. Geophys. Res. 93, 6463–6476 (1988).
(1983). 44. Bonadonna, C., Biass, S. & Costa, A. Physical characterization of explosive
16. Asano, G. Several discoveries obtained from the Hakutosan (Baegdusan). volcanic eruptions based on tephra deposits: propagation of uncertainties and
Expedition in 1942–1943. Mineral Geol. I, 228–231 (1948). sensitivity analysis. J. Volcanol. Geotherm. Res. 296, 80–100 (2015).
17. McLean, D. et al. Identification of the Changbaishan ‘Millennium’(B-Tm) 45. Constantinescu, R. et al. Uncertainty quantification of eruption source
eruption deposit in the Lake Suigetsu (SG06) sedimentary archive, Japan: parameters estimated from tephra fall deposits. Geophys. Res. Lett. 49,
synchronisation of hemispheric-wide palaeoclimate archives. Quat. Sci. Rev. e2021GL097425 (2022).
150, 301–307 (2016). 46. Bonasia, R., Macedonio, G., Costa, A., Mele, D. & Sulpizio, R. Numerical
18. Machida, H., Moriwaki, H. & Zhao, D. The recent major eruption of Changbai inversion and analysis of tephra fallout deposits from the 472 AD sub-Plinian
volcano and its environmental effects. Geogr. Rep. Tokyo Metrop. Univ. 25, eruption at vesuvius (Italy) through a new best-fit procedure. J. Volcanol.
1–20 (1990). Geotherm. Res. 189, 238–246 (2010).
19. Zhang, M. et al. The intraplate Changbaishan volcanic field (China/North 47. Volentik, A. C. M., Bonadonna, C., Connor, C. B., Connor, L. J. & Rosi, M.
Korea): a review on eruptive history, magma genesis, geodynamic significance, Modeling tephra dispersal in absence of wind: Insights from the climactic
recent dynamics and potential hazards. Earth Sci. Rev. 187, 19–52 (2018). phase of the 2450 BP Plinian eruption of Pululagua volcano (Ecuador). J.
20. Yun, S. H., Won, C. K. & Lee, M. W. Cenozoic volcanic activity and Volcanol. Geotherm. Res. 193, 117–136 (2010).
petrochemistry of volcanic rocks in the Mt. Paektu area. J. Geol. Soc. Korea 29, 48. Mingari, L. et al. Data assimilation of volcanic aerosol observations using
291–307 (1993). FALL3D + PDAF. Atmos. Chem. Phys. 22, 1773–1792 (2022).
21. Pan, B., de Silva, S. L., Xu, J., Liu, S. & Xu, D. Late Pleistocene to present day 49. Bonadonna, C. & Costa, A. Plume height, volume, and classification of
eruptive history of the Changbaishan-Tianchi Volcano China/DPRK: new volcanic explosive eruptions based on the Weibull function. Bull. Volcanol. 75,
field, geochronological and chemical constraints. J. Volcanol. Geotherm. Res. 1–19 (2013).
399, 106870 (2020). 50. Holasek, R., Self, S. & Woods, A. Satellite observations and interpretation of
22. Ikehara, K. Late Quaternary seasonal sea-ice history of the northern Japan Sea. the 1991 Mount Pinatubo eruption plumes. J. Geophys. Res. 101, 27,635–27,
J. Oceanog. 59, 585–593 (2003). 655 (1996).
23. Derkachev, A. N. et al. Tephra layers of large explosive eruptions of 51. Costa, A., Folch, A. & Macedonio, G. Density-driven transport in the umbrella
Baitoushan/Changbaishan Volcano in the Japan Sea sediments. Quat. Int. 519, region of volcanic clouds: Implications for tephra dispersion models. Geophys.
200–214 (2019). Res. Lett. 40, 4823–4827 (2013).
24. McLean, D. et al. Refining the eruptive history of Ulleungdo and 52. Scaillet, B. & Oppenheimer, C. Reassessment of the sulfur and halogen
Changbaishan volcanoes (East Asia) over the last 86 kyrs using distal emissions from the Millennium Eruption of Changbaishan (Paektu) volcano.
sedimentary records. J. Volcanol. Geotherm. Res. 389, 301–307 (2020). J. Volcanol. Geotherm. Res. 442, 107909 (2023).
25. Yun, S. H. Volcanological interpretation of historical eruptions of Mt. 53. Gerlach, T. M., Westrich, H. R. & Symonds, R. B. Preeruption vapor in
Baekdusan volcano. J. Korean Earth Sci. Soc. 34, 456–469 (2013). magma of the climactic Mount Pinatubo eruption: source of the giant
26. Wei, H. et al. Timescale and evolution of the intracontinental Tianchi volcanic stratospheric sulfur dioxide cloud. in Fire and Mud: Eruptions and Lahars of
shield and ignimbrite-forming eruption, Changbaishan, Northeast China. Mount Pinatubo, Philippines (eds Newhall, N. G. & Punongbayan, R. S.)
Lithos 96, 315–324 (2007). 415–433 (University of Washington Press, 1996).
27. Miyake, F., Nagaya, K., Masuda, K. & Nakamura, T. A signature of cosmic-ray 54. Sigl, M., Toohey, M., McConnell, J. R., Cole-Dai, J. & Severi, M. Volcanic
increase in AD 774–775 from tree rings in Japan. Nature 486, 240–242 stratospheric sulfur injections and aerosol optical depth during the Holocene
(2012). (past 11 500 years) from a bipolar ice-core array. Earth Syst. Sci. Data 14,
28. Sigl, M. et al. Timing and climate forcing of volcanic eruptions for the past 3167–3196 (2022).
2,500 years. Nature 523, 543–549 (2015). 55. Kravitz, B. & Robock, A. Climate effects of high‐latitude volcanic eruptions:
29. Chǒng Inji, K. Seoul National University Library [in Chinese, 高丽史] (1451). Role of the time of year. J. Geophys. Res. 116, D01105 (2011).
30. Chen, Z. Q. & Chen, Z. Identifying references to volcanic eruptions in Chinese 56. Staunton-Sykes, J. et al. Co-emission of volcanic sulfur and halogens amplifies
historical records. Geol. Soc. Lond. Spec. Publ. 510, 271–289 (2021). volcanic effective radiative forcing. Atmos. Chem. Phys. 21, 9009–9029 (2021).
31. Tsuboi, K. 坪井九馬三. Kohfukuji Chronicle. Handwritten copy from 57. Wade, D. C. et al. Reconciling the climate and ozone response to the 1257 CE
National Diet Library Collection. https://dl.ndl.go.jp/. (1908). Mount Samalas eruption. Proc. Natl Acad. Sci. 117, 26651–26659 (2020).
32. Fujiwara Tadahira [藤原忠平] 880–949, Teishin kō ki (Journal of Fujiwara 58. Prata, A. T., Mingari, L., Folch, A., Macedonio, G. & Costa, A. FALL3D-8.0: a
Tadahira) [貞信公記]. in Dai Nihon kokiroku [大日本古記録], Vol. 1 (eds computational model for atmospheric transport and deposition of particles,
Tokyo Daigaku Shiryō Hensanjo) (Iwanami Shoten, 1952–1986). Handwritten aerosols and radionuclides—part 2: model validation. Geosci. Model Dev. 14,
copy from National Diet Library Collection. https://dl.ndl.go.jp/. 409–436 (2021).
33. Horn, S. & Schmincke, H.-U. Volatile emission during the eruption of 59. Costa, A. et al. Quantifying volcanic ash dispersal and impact from
Baitoushan Volcano (China/North Korea) ca. 969 AD. Bull. Volcanol. 61, Campanian Ignimbrite super-eruption. Geophys. Res. Lett. 39, L10310 (2012).
537–555 (2000). 60. Macedonio, G., Costa, A. & Longo, A. A computer model for volcanic ash
34. Jin, B. & Zhang, X. Y. Researching Volcanic Geology in Changbai Mt. (in fallout and assessment of subsequent hazard. Comput. Geosci. 31, 837–845
Chinese) (Northeast Korean Nationality Education Press, 1994). (2005).

12 COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:6 | https://doi.org/10.1038/s43247-023-01162-0 | www.nature.com/commsenv

COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-01162-0 ARTICLE

61. Pfeiffer, T., Costa, A. & Macedonio, G. A model for the numerical simulation Author contributions
of tephra fall deposits. J. Volcanol. Geotherm. Res. 140, 237–294 (2005). Conceptualisation: A.C. and V.C.S.; Methodology: A.C., L.M., and V.C.S.; Software: L.M.,
62. Costa, A., Pioli, L. & Bonadonna, C. Assessing tephra total grain-size A.F., G.M., and A.C.; Resources: L.M., V.C.S., D.M., J.L., and S.H.Y.; Writing—original
distribution: insights from field data analysis. Earth Planet. Sci. Lett. 443, draft: A.C., L.M., and V.C.S.; Writing—review and editing: A.C., L.M., V.C.S., G.M.,
90–107 (2016). D.M., A.F., J.L., and S.H.Y.; Visualisation: L.M.; Supervision: A.C. and S.H.Y.; Funding
63. Marti, A., Folch, A., Costa, A. & Engwell, S. Reconstructing the plinian and acquisition: A.C. and A.F. All authors have read and approved the final version of the
coignimbrite sources of large volcanic eruptions: a novel approach for the manuscript.
Campanian Ignimbrite. Sci. Rep. 6, 21220 (2016).
64. Yi, J. et al. Modeling the multi-level plumbing system of the Changbaishan
caldera from geochemical, mineralogical, Sr-Nd isotopic and integrated
Competing interests
The authors declare no competing interests.
geophysical data. Geosci. Front. 12, 101171 (2021).
65. Bonadonna, C. & Phillips, J. C. Sedimentation from strong volcanic plumes. J.
Geophys. Res. 108, 2340 (2003). Additional information
66. Mele, D. et al. Total grain size distribution of components of fallout deposits Supplementary information The online version contains supplementary material
and implications for magma fragmentation mechanisms: examples from available at https://doi.org/10.1038/s43247-023-01162-0.
Campi Flegrei caldera (Italy). Bull. Volcanol. 82, 31 (2020).
67. Costa, A., Folch, A. & Macedonio, G. A model for wet aggregation of ash Correspondence and requests for materials should be addressed to Antonio Costa.
particles in volcanic plumes and clouds: I. Theoretical formulation. J. Geophys.
Res. 115, B09201 (2010). Peer review information Communications Earth & Environment thanks Céline Vidal
68. Cornell, W., Carey, S. & Sigurdsson, H. Computer simulation and transport of and the other, anonymous, reviewer(s) for their contribution to the peer review of this
the Campanian Y−5 ash. J. Volcanol. Geotherm. Res. 17, 89–109 (1983). work. Primary Handling Editors: Joe Aslin and Carolina Ortiz Guerrero. A peer review
69. McKay, M. D., Beckman, R. J. & Conover, W. J. A comparison of three file is available
methods for selecting values of input variables in the analysis of output from a
computer code. Technometrics 21, 239–245 (1979). Reprints and permission information is available at http://www.nature.com/reprints
70. Suzuki, T. A. theoretical model for dispersion of tephra. In Arc Volcanism:
Physics and Tectonics, (eds Shimozuru, D. & Yokoyama, I.) pp 93–113 (Terra Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in
Scientific Publishing Company (TERRAPUB), Tokyo, 1983). published maps and institutional affiliations.
71. Degruyter, W. & Bonadonna, C. Improving on mass flow rate estimates of
volcanic eruptions. Geophys. Res. Lett. 39, L16308 (2012).
Open Access This article is licensed under a Creative Commons
Attribution 4.0 International License, which permits use, sharing,
Acknowledgements adaptation, distribution and reproduction in any medium or format, as long as you give
A.C., G.M. and L.M. thank the Istituto Nazionale di Geofisica e Vulcanologia, Italy, grant appropriate credit to the original author(s) and the source, provide a link to the Creative
Progetto INGV Pianeta Dinamico (code CUP D53J19000170001) funded by Italian Commons license, and indicate if changes were made. The images or other third party
Ministry MIUR (“Fondo Finalizzato al rilancio degli investimenti delle amministrazioni material in this article are included in the article’s Creative Commons license, unless
centrali dello Stato e allo sviluppo del Paese”, legge 145/2018). This work has been indicated otherwise in a credit line to the material. If material is not included in the
partially funded by the Centre of Excellence for Exascale in Solid Earth (ChEESE; article’s Creative Commons license and your intended use is not permitted by statutory
EuroHPC Grant Agreement No. 101093038), the project “A Digital Twin for GEO- regulation or exceeds the permitted use, you will need to obtain permission directly from
physical extremes” (DT-GEO; Horizon Europe Grant Agreement No. 101058129), and the copyright holder. To view a copy of this license, visit http://creativecommons.org/
the Meteorological/Earthquake See-At Technology Development Research Grant licenses/by/4.0/.
KMI2018-02710. We acknowledge European Centre for Medium-range Weather Fore-
cast (ECMWF) for furnishing the meteorological re-analysis ERA-5 data (no permissions
were required). We thank two anonymous reviewers and Céline Vidal for constructive © The Author(s) 2024
criticism and feedback that improved this manuscript.

COMMUNICATIONS EARTH & ENVIRONMENT | (2024)5:6 | https://doi.org/10.1038/s43247-023-01162-0 | www.nature.com/commsenv 13