Assessing Volcanic Hazard and Exposure at Obscure Volcanic Fields: A Case Study from the Bolaven Volcanic Field, Laos

Southeast Asia hosts a large number of active and well-studied volcanoes, the majority of which are 14 located in Indonesia and the Philippines. Northern Southeast Asia (Myanmar, Cambodia, Laos, 15 Thailand and Vietnam) also hosts volcanoes that for several reasons (post-World War II conflicts, poor 16 accessibility due to dense vegetation, no known historical activity) have been little studied. Systematic 17 assessments of the threat these volcanoes pose to resident populations do not exist, despite evidence 18 of numerous eruptions through the late Pleistocene and likely even during the Holocene. A recent 19 study that inferred the location of the Australasian meteorite impact (which produced the largest 20 known tektite strewn field on Earth) beneath the Bolaven Volcanic Field in southern Laos provided a 21 wealth of data for that volcanic field, in particular, mapping of vents and flows, and their absolute 22 ages. Building upon this foundation, we used the Bolaven Volcanic Field as a case study for assessing 23 the potential exposure of populations and infrastructure to lava flows during future eruptions there. 24 Our study uses remote sensing of past flows, lava-flow simulations and open-access exposure data, to 25 assess hazards and exposure. Our results show that future vents are most likely to occur in a N-S band 26 atop the Bolaven Plateau, with some flows channelled into canyons that spill down the plateau flanks 27 onto lower plains that support more populated areas such as the provincial centre, Pakse. Our Black Point – Pyroclasts of a Surtseyan show no change during edifice growth to the from


Introduction
historically active volcanoes, including Mount Popa, Myanmar, which last erupted in 442 BCE). Both 48 works highlighted the presence of young lava flows and well preserved scoria cones and crater lakes. 49 The young (and even historical) ages of large lava flows encouraged us to conduct a more-focused 50 study on the hazards represented by, exposures to, and risks from these volcanoes. 51 There are several reasons why such focused studies have not taken place previously. These include 52 1) the dangers of field work during and subsequent to the conflicts of the second half of the 20 th 53 century, including the presence of still-unexploded ordnance; 2) the limited access due to dense 54 tropical vegetation, as well as poor exposures; and 3) the scant historical eruptions, compared to the 55 abundant historical eruptions along the volcanic arcs of Indonesia and The Philippines. Thus, the 56 potential for future eruptions of these volcanic fields to affect growing and thriving populations, 57 agriculture and tourism is under-appreciated and unknown. 58 For this study, we take advantage of a recent and significant increase in knowledge about one of 59 for lava flows from 37 locations. We utilize these data and additional field data and 30 additional 68 geochronological dates to provide this first volcanic hazard and exposure assessment for the region.  Jiménez et al., 2020). In most efforts, geoscientists have conducted volcanic hazard assessments on well-studied and easily accessible volcanoes. The BVF is neither well-studied nor easily accessible. 82 Here we utilize a combination of field, remote-sensing, and numerical-modelling techniques that can 83 be applied to other understudied, volcanic fields with difficult access. 84 The primary hazard associated with the BVF is lava flows. Previous studies that considered 85 exposure to lava flows (e.g. Bonne  The BVF has no historical eruptive record. It is clear, however, that there are key assets that could 90 be exposed to future lava flows. The potential infrastructure damage and economic loss from an 91 eruption of the BVF may be significant. A quantitative assessment of exposure to the BVF could help 92 Laotian communities and governments and their partners in planning for future volcanic crises. The 93 national government recently requested such an assessment as part of a geohazards collaboration 94 between the government of Laos and the Earth Observatory of Singapore, through the CCOP (the SE 95 map in Fig.  2 and detailed map in a data repository at 181 https://researchdata.ntu.edu.sg/privateurl.xhtml?token=605d1696-0a8d-47cf-afeb-5e18398a6ef2). 182 The updated map includes new field observations, 30 new 40 Ar-39 Ar dates of lava flows, and new 183 interpretations based on digital imagery from ©Google Earth Pro and a geospatial analysis using Esri® 184 ArcMap 10.7.1 (Fig. 3). Field observations included: 1) The use of vegetation as a proxy to determine 185 the relative age of the very youngest flows (i.e. an immature virgin forest implies a very young age), 186 before performing 40 Ar-39 Ar dating; 2) Geomorphological evidence (e.g. differential erosion of 187 neighbouring flows indicating different relative age, or evidence of rivers displaced by lava flows); and 188 3) Different macro-scale flow textures (e.g. colour, relative abundance of minerals and mineral type, 189 level of weathering). Satellite imagery was mostly used in combination with field observations. 190 ArcMap 10.7.1 was used to create topographic contours from a Digital Elevation Model (Shuttle Radar  191 Topography Mission, 30-m resolution), in order to identify, where possible, different flows (even 192 within the same flow complex) and topographic features associated with them (e.g. scoria cones). 193 Scoria cones were mapped based on (i) height (at least a few tens of meters high), (ii) morphology, 194 which also provided clues on their relative age (i.e. conical, horseshoe-like, flat-topped, with the first 195 and the last being relatively younger and older respectively), and (iii) location respect to the nearby 196 lava flow(s) (surrounded or incorporated by/in the nearby flow(s), to define if they were older or the same age of that flow). Besides scoria cones, other positive-relief topographical features were present, 198 and were mapped as "mounds". A common characteristic of these mounds is the relatively low H/W 199 ratio (height less than 10 m, occasionally up to 30 m; width from a few hundred meters up to a few 200 kilometres Elements of a KDE are the kernel function and the bandwidth. The kernel function defines the 224 probability of future vent formation at locations within a certain region, and can be of different 225 statistical types (e.g. Gaussian, Cauchy, Epanechnikov, Triangular, Uniform, Triweight, Quartic); the 226 bandwidth is the search radius within which the density is calculated. The latter, in particular, is more 227 sensitive to the resulting output relative to the type of kernel function (Martin et al., 2004), and 228 therefore requires careful selection. We performed the KDE with ArcMap 10.7.1, which uses a Quartic 229 kernel function (Silverman, 1998) and provides a default bandwidth (defined as search radius in the 230 software) of ~7.4 km, computed specifically for the input dataset (distribution and number of scoria 231 cones across the BVF) using a spatial variant of Silverman's Rule of Thumb (Silverman, 1998). This 232 method, compared to more classic spatial density estimations, is weighted for spatial outliers, which 233 pertains to some vents on the BVF (Fig. 4). 234 As input data for the KDE on the BVF, we selected vents based upon two main criteria: 1) Our 235 confidence that the source is indeed a vent; and 2) The vent age is younger or equal to 790 ka (the age 236 of the Australasian impact crater). For the first criterion, we binned eruptive vents on the BVF into 237 three types: Scoria cones, mounds and fissures. Only scoria cones (n= 76) were considered for the 238 estimation of the vent spatial distribution on the BVF. We ignored the many low-relief mounds and 239 fissures, because we know from ground-truthing some of these that the dense and tall vegetation 240 common in the region is the cause of some of these features in the DEM. Moreover, the relative 241 vertical accuracy of the SRTM is around 6-m at these latitudes (Brown et al., 2005), about the same as 242 the height of these features. We constrained vents to those younger than the impact crater age as a 243 reference for the volcanic hazard assessment of the BVF for two main reasons: i) The Australasian 244 meteorite impact represented an important event in the volcanic history of the Bolaven Plateau. 245 Although volcanism began long before the impact, there are hints that the event may have affected 246 the rate of melt production (Sieh et al., 2019), as has been inferred for other large impacts (Jones, 247 2005), and as changes in the geochemical composition of the basalts suggests (Herrin et al., in prep.).

256
The output of the KDE is a vent spatial-density map with seven classes of visualization, each one 257 corresponding to a range within the density field, expressed as the number of scoria cones per unit 258 area. We selected 30-m as the cell size for the output map, matching the resolution of the SRTM. 259

Defining eruption source parameters 276
The range of eruption source parameters used for our lava flow simulations are listed in Table 1,  277 and include thickness, volume (total volume to be erupted) and pulse volume (maximum volume for 278 each pulse). In order to choose adequate parameters, we used a combination of field data from the 279 BVF (where measurable or available) and from analogue volcanic fields. We lack a complete record of  Pulse volume (m 3 ) 6.80×10 4 4.59×10 5 n/a n/a 291

Model inputs 292
Other inputs for MOLASSES included a Digital Elevation Model (DEM) of the area and previous vent 293 locations. We used the Shuttle Radar Topography Mission (SRTM) DEM, which is available at 1 arc-294 second resolution (30m × 30m at BVF), and which we cropped to a total extent of ~12,000 km 2 (102 295 km × 118 km grid) to cover the BVF. The locations of the new vents were stochastically sampled 296 according to the vent spatial density map (Fig. 4). One lava flow was simulated for each vent, with the 297 eruption volume and average lava flow thickness stochastically sampled from within the range 298 considered (Table 1) Using resolutions finer than 30-m for the BVF may provide better results, however, two main issues 320 would be the availability of such DEMs and the computation time required to run 10,000 simulations. 321 Therefore, based on these considerations we decided to proceed with the 30-m resolution DEM. 322

Model outputs 323
MOLASSES outputs a hit intensity map, on which the intensity is the number of times that each grid 324 cell is inundated (hit) by lava flows. By dividing the number of hits in each grid cell by the total number 325 of simulations performed, we obtained a conditional probability of inundation for each grid cell. A 326 conditional probability assumes that an eruption has occurred somewhere in the BVF. 327 We defined three hazard zones, coloured yellow, orange and red, based on the probability of 328 inundation, at the 90 th (1-78 hits), 50 th (78-390 hits) and 10 th (390-779 hits) percentiles. This choice 329 was somewhat arbitrary but taken to reflect the distributions of hits across the BVF, into areas with 330 relatively low, medium and high probability of inundation. 331 332 4.6. Population, infrastructure and land cover exposure 333 We considered population, power lines, power stations, dams, roads, and land cover (i.e. forested 334 and vegetated areas, croplands, built-up areas and water bodies) as critical elements to consider for 335 the BVF because of the intense use of land for agriculture and hydropower production and supply. 336 Other aspects of land use, such as building type and purpose, were not considered here due to the 337 lack and/or reliability of data , or because the exposure to future lava flows was considered very low 338 (e.g. Pakse international Airport). . Moreover, the region of the BVF is predominantly rural; these characteristics make population 346 estimates less reliable when geospatial elements such as nightlights are used as an indication of 347 population density and size (Small et al., 2005). Therefore, we decided to refer to three different, 348 widely used, free and relatively up-to-date population datasets, with different spatial resolutions, in 349 order to provide an indication of the uncertainty in our estimate of people exposed to a future 350 eruption. We used WorldPop2020   Table 2). The former consists of roads for public use (e.g. country's 375 roads, links between cities), the latter includes roads for private or semi-private use (e.g. access to 376 housing, industrial or agricultural use). Note that for the purpose of this paper, roads for exclusive 377 pedestrian use were not included.

Secondary link
Used to identify slip roads/ramps and "channelised" at-grade turning lanes connecting the through carriageways/through lanes of a Secondary other minor roadways

Tertiary
Roads connecting smaller settlements, and within large settlements for roads connecting local centers. In terms of the transportation network, OpenStreetMap "tertiary" roads commonly also connect minor streets to more major roads.

Unclassified
Minor public roads typically at the lowest level of the interconnecting grid network

Residential
Roads for accessing residential areas and in residential areas but not normally used as through routes

Service
Roads for access to a building, service station, beach, campsite, industrial estate, business park, etc.

Track
Roads mostly for agricultural use, forest tracks etc. Often unpaved (unsealed) but may be paved tracks Suitable for two-track vehicles, such as tractors or jeeps a OSM= Open Street Map 391 392

Lava flows 398
Post-meteorite impact lava flows vary in length between a few hundred meters and a few tens of 399 kilometres (up to ~50 km, northern flow complex), for a total areal extent of about 3900 km 2 . basalts on the field (Fig. 3). 409

Scoria cones 410
Among the 67 available dates from the BVF, 47 are from post-meteorite impact lavas ( Fig. 2 and  411 Table S1, supplementary material). Two of these are from lava within scoria cones. Geochronological 412 and geochemical data from scoria cones are sparse, because they are commonly in thick, untracked 413 jungle and their clastic deposits are more highly weathered than lava flows (only one scoria cone was 414 successfully dated, ~1.5 km E of Paksong, 200 ± 9 ka). Therefore, most scoria cones were either 415 assigned an absolute age, based on 40 Ar-39 Ar dates of lava flows that erupted from that cone (e.g. Fig.  416 3e), or assigned an age interval (if an absolute age was not available for any of the flows that emanated 417 from beneath that cone). Out of the 97 scoria cones mapped (Fig. 2), 76 were inferred to be younger 418 than 790 ka, whereas 21 were either older than 790 ka or indeterminate in age/age-interval. 419 Here we assume that each vent was produced by a single eruptive event (following the definition 420 of vents vs events as in Gallant et al., 2018). This assumption was based on the lack of field information 421 for most scoria cones (lack of fresh scoria cones deposits and/or lack of accessibility), resulting in 422 missing information usually needed to link an eruptive source to its deposits; for example, use of 423 geochemical compositions and absolute ages of scoria cones, to match (or not match) them with 424 nearby lava flows and other scoria cones. More discussion around the issue of vents vs events is 425 reported in Section 6. 426 427

Vent spatial density and probabilistic inundation map of the BVF 428
The approximately 30% of scoria cones within the highest density fields (5 th to 33 th percentile) form 429 a N-S band across the NW portion of the plateau (Fig. 4). Some isolated clusters of lower density are 430 to the NE (~15 km south of Salavan), to the W (~10km north-east of Pakse), and about 35 km N of 431 Pakse. Curiously, the area with the highest concentration of vents is nearly coincident with the inferred 432 location of the Australasian impact crater (Fig. 4). using LandScan2019, are 80% in the yellow zone, 16% in the orange zone, and 4% in the red zone 450 (Table 3). 451

Infrastructure 452
There are 2678 km of roads within the three potentially inundated areas (Fig. 5). Of these, 1479 453 km (55%) are within the yellow zone, 950 km (35%) within the orange zone, and 250 km (10%) are 454 within the red zone. In terms of road type, the balance is slightly in favour of Type-1 roads for yellow 455 and orange zones (55% and 59% respectively), and roughly the same for the red zone (Type-1= 49%, 456 Type-2= 51%). This reflects the underlying proportion of roads, which are approximately 56% Type-1 457 and 44% Type-2 across our study area (see also Fig. 5). 458 Across the entire hazard area are 416 km of power lines. Of these, 259 km (62%) are in the yellow 459 zone, 120 km in the orange zone (29%), and 37 km (9%) in the red zone. Along these power lines, there 460 are two power stations (both in the yellow zone) and one dam (orange zone). Another dam is located 461 in the orange zone (

Land Cover 477
Across all hazard zones (Fig. 6), 4996 km 2 (89%) of land is covered in forest and sparse vegetation. 478 Cropland covers 514 km 2 (9%), built-up areas cover 73 km 2 (1%), and permanent and seasonal water 479 cover 15 km 2 (<1%). In particular, the yellow zone encompasses 3060 km 2 (87%) of forest and sparse 480 vegetation, 387 km 2 (11%) of cropland, 48 km 2 (1%) of built-up areas, and 13 km 2 (>1%) of water. The Volcanic risk assessments require information on several aspects, such as volcanological setting 492 (e.g. central volcano, volcanic field), eruptive style (e.g. explosive, effusive), hazard type (e.g. ash 493 fallout, pyroclastic flows, lava fountains, lava flows, lahars), and the geographic context of the area 494 assessed (how the hazards can impact the local population and/or their activities). Their robustness 495 often depends strongly upon the availability of data for that particular area (reliable and up-to-date 496 geological, hazard, exposure and vulnerability information). Here we did not conduct any vulnerability 497 assessment due to the lack of information. Below, we discuss the results and limitations of the present 498 research, based on evaluating likely hazard, exposure and potential impacts in the BVF area in case of 499 renewed volcanism, the first such assessment in northern SE Asia for a volcanic field. 500 Understanding the distribution of eruptive centres within a volcanic field may help understand the 506 likely location of future activity at that field, hence allowing more accurate hazard and exposure 507 assessments. The location of the post-impact scoria cones on the BVF, hence the area with the highest 508 probability of future vent opening, seems related to the impact crater structure rather than to the 509 local tectonic stress. 510 From a tectonic perspective, a report from the Ministry of Energy and Mines (Japan International 511

Cooperation
Agency, 2008) and mapping effort from this work 512 (https://researchdata.ntu.edu.sg/privateurl.xhtml?token=605d1696-0a8d-47cf-afeb-5e18398a6ef2), 513 highlighted the presence of several tectonic structures (a syncline crossing the plateau, a 514 syncline/anticline pair just NE of the Plateau, and several reverse faults to the E of the plateau), whose 515 σ1 is oriented NE-SW. Field evidence and Ar-Ar dating of a nearby lava flow indicate that these 516 structures (particularly the syncline/anticline pair NE of the Plateau) are younger than 2.9 Ma, and 517 they likely represent the latest deformation pattern in this region. If there was a direct link between 518 this local stress field and the pattern of scoria cones, then it would be parallel to the σ1 (Nakamura, 519 1977), which is not the case here, given the N-S orientation of the scoria cones. 520 On the other hand, the N-S aligned post-impact scoria cones on the BVF, besides being 521 encompassed within the inferred location of the Australasian impact crater, have an inferred age 522 lower than the impact itself, with pre-impact scoria cones showing a more widespread distribution 523 across the field (Fig. 4). This may suggest an existing link between the meteorite impact and the post-524 impact volcanism, as hypothesized for other large-scale meteorite impacts (Jones, 2005), or in part 525 may reflect a bias from the lack of exposure of older scoria cones (buried or eroded), potentially hiding 526 older (or unknown) regimes. 527 Although additional data are needed (including a complete dataset of seismic tomography, 528 geochemical data, and absolute ages of scoria cones) to better constrain the reasoning behind this distribution, our scoria cones spatial distribution analysis provides a first-time indication of the 530 location of future volcanic activity on the BVF, which largely coincides with the location of the largest 531 known young meteorite impact on Earth. 532 533

Exposure analysis 534
In this section we discuss exposure and impacts on populations, infrastructures and landcover, 535 based on our lava flow simulations and mapping. 536 The population exposure results show that there is some variation among all datasets. In particular, 537 the highest variation is observed for LandScan, where the total count is between 13% and 24% smaller 538 than that for WorldPop and GHS-POP (Table 3) respectively, and also different for individual hazard 539 zones (20% higher for yellow zone, 17% smaller for orange zone, and 3% smaller for red zone). GHS-540 POP has the highest total population count (~358,000), but if we consider the fractions for each hazard 541 zone, WorldPop and GHS-POP share the same % of people. All three datasets are based on census data 542 and geospatial information. Decadal censuses are conducted in Laos, with the last conducted in 2015, 543 hence relatively recent. A major difference among the considered datasets is the spatial resolution, 544 with LandScan having the coarsest one. In addition, as anticipated in Section 4, LandScan has been 545 found to underestimate population counts at the transition between urban and rural areas (Calka and 546 Bielecka (2019) and the lower LandScan population count on the BVF may reflect the presence of 547 sparse villages in the region, surrounded by forests and croplands (Fig. 6). WorldPop and GHS-POP 548 may present some limitations as well, however, the results are relatively consistent when total counts 549 and fractions across the different hazard zones are compared. For the purpose of this study, the 550 application of largely used, and freely available population datasets allowed us to detect the minimum 551 number of people potentially exposed to lava flows in case of an eruption on the BVF. impact based on factors such as the characteristics of the flow (e.g. viscosity, advancement rate, 561 thickness), topography and distance between the exposed assets and the vent(s), people's behaviour (e.g. maintaining a safe distance from the flows), quality and type of infrastructure, and response 563 strategies (e.g. redirection of lava flows to less exposed/vulnerable areas). Education of the population 564 to the likely hazards posed by lava flows, and their long-term impacts can help communities adapt and 565 recover . 566 Cities on and around the BVF have calculated conditional probability of inundation that vary from 567 0.01% to 7.79% across the field (Fig. 5). Pakse is the largest and most populated city in southern Laos, 568 and 40 Ar-39 Ar dating shows that the majority of the city is built upon a lava flow ~180 ka old, and a flow 569 likely <40 ka is nearby. This low spatial density of past flows puts Pakse within the least exposed zone, 570 our Yellow Zone (Fig. 5). It has about a 0.05% likelihood of being reached by lava flows during any 571 future eruption. That is to say, there is one chance in 2,000 that any future flow will enter the city. 572 Paksong, the second largest city in the area, sits within the Red Zone. Although the flow that it sits 573 upon is about the same age as the one underlying Pakse (~185 ka), Paksong is within a few kilometers 574 of many young flows that range in age from ~75 to ~200 ka. The likelihood of a future flow reaching 575 Paksong is ~6% (or about one chance in 17), about a hundred times higher than one reaching Pakse. 576 The city of Salavan is north of the BVF and is built on young fluvial sediments rather than lava flows. 577 Our simulations thus show a conditional probability of inundation of 0%. That is not to say that there 578 is absolutely no possibility of inundation by lava, but the likelihood is vanishingly small. Salavan is more 579 likely to be affected by other aspects of an eruption. For example, a flow down the north side of the 580 plateau could reach and dam the Xedon River, flooding parts of the plain around and including the 581 city. 582 People can be significantly affected by damage to infrastructure. Disruption of power production 583 and distribution could result in loss of access to power for most areas in the Champasak, Salavan, and 584 Xekong provinces. Here we only consider hydropower, for the reasons explained in Section 2, which 585 is produced by the operational dams in the area, and distributed through a network of power lines 586 (Fig. 5). Eruption-related damage to one of the two dams located in the orange zone, which is directly 587 connected to one of the power lines, may create a disruption in the immediate power production and 588 supply. Power lines on and around the BVF cover the entire width of the potentially inundated area, 589 therefore, regardless of where the damage will take place, it may affect other areas of the plateau, 590 potentially also affecting urban sites located outside of this area (e.g. Salavan). In addition, most of 591 the power produced in Laos is exported to Thailand; therefore, international trading can undergo an 592 impact as well. Another key infrastructure that can suffer an impact, affecting people, is roads. This 593 can happen through restricted or blocked access for: 1) delivery of essential goods (e.g. food, 594 medicines) to the areas affected; 2) emergency response or maintenance vehicles, for example to 595 tackle wildfires that may be initiated by the flows or to repair damage to power infrastructure; 3) people, who may not be able to reach their work sites or visit family members in nearby cities or 597 villages who need assistance; 4) trading with nearby countries (Thailand, Cambodia and Vietnam). 598 Another key asset to consider in exposure assessments is land cover. Our analysis of the potentially 599 inundated area shows an obvious predominance of forest and other types of natural vegetation, 600 followed by croplands, built-up areas and water bodies (Fig. 6). An eruption on the Bolaven plateau 601 affecting them can result in impacts for the local and national economy. Forests, for example, may be 602 ignited by lava flows (e.g. Ainsworth and Kauffman, 2009;Harris, 2015), particularly in the drier 603 season, and fires can propagate over a large area, potentially damaging infrastructures and natural 604 historical sites. Cropland and built-up land classes, which seem to be linked across the whole field 605 (built-up areas being surrounded or adjacent to areas with high cropland fraction) double from orange 606 zone to red zone (5% to 10% for cropland, and 1% to 2% for built-up areas), despite the red hazard 607 zone being ~6 times smaller than the orange hazard zone; this suggests that the high area of the 608 plateau (also the red zone: Fig 6) has the most favourable conditions for agricultural use. Although 609 cropland was not subdivided into different agricultural classes here (lack of reliable and/or up-to-date 610 data), coffee production on the Bolaven plateau represents about 95% of the total amount produced 611 in Laos (Toro, 2012). A future eruption on the Bolaven plateau could therefore affect the socio-612 economic wellbeing of those in the region, who rely on coffee as their source of income. This in turn 613 can impact the economy of the whole country (Toro, 2012), either directly, if coffee plantations are 614 inundated by lava, or indirectly, if roads linking cities or countries are inundated by lava and 615 inaccessible for months/years affecting transport to major cities for coffee processing and export. 616 Official GDP data from Laos are available for the agriculture sector as a whole, but not available for 617 individual sub-sectors such as coffee. 618 In order to further evaluate the knock-on consequences from an effusive eruption on the BVF, up-619 to-date information about infrastructure and land cover is needed. 620 621 6.3. The role of external water in potential explosive activity on the BVF 622 Explosive interaction between magma and water (here broadly referred to as phreatomagmatic 623 activity), is a potentially hazardous volcanic phenomenon known to occur in a wide range of 624 environmental settings (e.g. Thorarinsson, 1964 found no direct evidence of past explosive magma-water interaction on the BVF, the abundance of 628 water in this area may lead to such activity in case of eruption. Below we address three different