Assessing volcanic hazard and exposure to lava flows at remote volcanic fields: a case study from the Bolaven Volcanic Field, Laos

Abstract

Southeast Asia is home to a large number of active and well-studied volcanoes, the majority of which are located in Indonesia and the Philippines. Northern Southeast Asia (Myanmar, Cambodia, Laos, Thailand and Vietnam) also hosts volcanoes that for several reasons (post-World War II conflicts, poor accessibility due to dense vegetation, no known historical activity) have been poorly studied. Systematic assessments of the threat these volcanoes pose to resident populations do not exist, despite evidence of numerous eruptions through the late Pleistocene and likely even during the Holocene. A recent study inferred the location of the Australasian meteorite impact to be beneath the Bolaven Volcanic Field in southern Laos; this study provided a wealth of data for the field: in particular, mapping of vents and flows, and their relative or absolute ages. The Bolaven Volcanic Field (16 Ma—< 40 ka) has a surface area of about 5000 km², contains nearly 100 scoria cones and more than 100 individual lava flows. Some lava flow systems are as long as 50 km, with thickness ranging from a few meters at the flow edges, up to > 50 m in some locations. Building upon this foundation, we used the Bolaven Volcanic Field as a case study for assessing the potential exposure of populations and infrastructure to lava flows during future effusive eruptions. Our study uses remote sensing to map past flows and vents (i.e. scoria cones), lava-flow simulations from new simulated vents, and open-access exposure data, to assess hazards and exposure. Our results show that future vents are most likely to occur in a N-S band atop the Bolaven plateau, with some flows channelling into canyons and spilling down the plateau flanks onto lower plains that support more populated areas such as the provincial centre, Pakse. Our exposure assessment suggests that around 300,000 people could experience socio-economic impacts from future lava flow inundations. The largest impacts would be on two of the main economic sectors in the region, agriculture and hydropower. The potential also exists for life-threatening explosions from interactions between magma and surface waters, which are abundant in the region. We estimate an average recurrence interval of approximately 10,400 years, based on information from lava flows and scoria cones.

Introduction

Southern Vietnam, northern Cambodia, eastern Thailand and southern Laos host many Late-Cenozoic monogenetic basaltic volcanic fields, scattered over about 50,000 km² (Fig. 1). Very few studies of these fields exist, and mainly focus on regional tectonics and possible magmatic sources. Barr and Macdonald (1981) argued that basalts from these volcanic fields comprise a large continental basaltic province, which includes the submarine volcano Ile des Cendres, ~ 120 km south-east of Vietnam (Figure S1 – supplementary material), which formed two ephemeral islands during its last eruption in 1923. They note, however, that the temporal and geochemical variability of the region’s basalts do not show any particular spatial pattern. The eruptive products range from alkaline to tholeiitic to olivine basalts, and span ages from about 0.5 Ma to 12 Ma (Barr and Macdonald 1981). They hypothesize that this unpatterned variability reflects the presence of discrete magma sources interacting at various times, contributing magmas from various degrees of partial melting in a region of complex tectonic interactions. Whitford-Stark (1987) conducted a similar study to characterize the basalts across a much larger swath of mainland Asia and including this region. They characterized the distribution, age, geochemistry and origin of the volcanics. These authors also list a few potentially historically active volcanoes, including Mount Popa, Myanmar (Figure S1), which was inferred to have erupted in 442 BCE (Whitford-Stark 1987). Both works highlighted the presence of young lava flows and well preserved scoria cones and crater lakes.

Fig. 1

Late Cenozoic volcanic fields in northern Southeast Asia. The Bolaven Volcanic Field is the northernmost of these fields, in southern Laos. The inset shows the Bolaven plateau (elevation information are from the Shuttle Radar Topography Mission – SRTM) and the Laotian provinces within which it resides

There are several reasons why such focused studies have not taken place previously. These include 1) the dangers of field work during and subsequent to the conflicts of the second half of the twentieth century, including the presence of still-unexploded ordnance; 2) the limited access due to dense tropical vegetation, as well as poor exposures; and 3) the scant historical eruptions, compared to the abundant historical eruptions along the volcanic arcs of Indonesia and the Philippines. Thus, the potential for future eruptions at these volcanic fields is under-appreciated and unknown, but given the large number of people living in northern Southeast Asia (over 100 million), and its worldwide economic relevance (e.g. export of coffee, tropical fruits, rubber, etc.), a quantitative exposure assessment to volcanic hazards is needed.

For this study, we take advantage of a recent and significant increase in knowledge about one of the volcanic fields in northern SE Asia, the Bolaven Volcanic Field (BVF) in southern Laos (Fig. 1). Lava flows of the volcanic field mostly occupy the broad summit area of the Bolaven plateau, and spill down its flanks on the western and northern parts, covering an area of about 5000 km². The field is underlain by Mesozoic fluvial/lacustrine mudstones and sandstones, which reach a maximum thickness of about 900 m on the plateau edges (exposed on the western, southern and eastern edges). A 7-year project on the Bolaven plateau combined geological, petrological, geochemical and geophysical studies to identify a 17 km by 13 km ~ 790 ka impact crater related to the tektites of the Australasian strewn field (Sieh et al. 2019). The crater lies wholly buried beneath the summit of the BVF, covered by post-impact lavas. In their quest to identify the crater location, Sieh et al. constructed a detailed geological map of the BVF, which included lava flows and scoria cones, and determined ⁴⁰Ar-³⁹Ar dates for lava flows from 37 locations. We utilize these data and additional field data and 30 new geochronological dates to provide this first volcanic hazard and exposure assessment for the region. This combined effort resulted in mapping of 76 scoria cones and over 100 individual lava flows, among other geological features in the area (full geological map available at https://researchdata.ntu.edu.sg/dataset.xhtml?persistentId=doi:10.21979/N9/HQDXRQ&faces-redirect=true), and there seems to be a link between the spatial distribution of the scoria cones and the meteorite impact site.

Volcanic hazard assessments take into account spatial and/or temporal information from past activity of a volcano or volcanic field to forecast likely future activity (e.g. Cappello et al. 2012; Connor et al. 2012; Gallant et al. 2018). By coupling information on the hazard with that on exposure (e.g. the number and distribution of people and infrastructure within the affected area), we can make an estimate of the range of potential consequences from a future eruption (Barsotti et al. 2018; Jiménez et al. 2019; Jenkins et al. 2022). Such estimates support local authorities in characterising and preparing for future volcanic crises (e.g. Orsi et al. 2004; Felpeto et al. 2007; Bevilacqua et al. 2015; Jiménez et al. 2020). In most efforts, geoscientists have conducted volcanic hazard assessments on well-studied and easily accessible volcanoes (Connor et al. 2000; Orsi et al. 2004; Hayes et al. 2018; Hopkins et al. 2021). A few studies have been carried out on less-accessible volcanic fields (e.g. Walker et al. 2011; Ureta et al. 2021; Sieron et al. 2021). The BVF is neither well-studied nor easily accessible. Here we utilise a combination of field, remote-sensing, and numerical-modelling techniques that can be applied to other understudied volcanic fields with access restrictions.

The primary hazard associated with the BVF is lava flows, given the lack of evidence for other possible hazards (e.g. presence of tephra). From a broader perspective, previous studies assessing exposure to lava flows have mainly focused on currently active volcanoes (e.g. Bonne et al. 2008; Bisson et al. 2009; Del Negro et al. 2020). While, only a few focused on historically dormant volcanoes (e.g. Deligne et al. 2017; Hayes et al. 2021). The BVF has no historical eruptive record. It is clear, however, that there are key assets that could be exposed to future lava flows, such as coffee plantations (a major businesses in Laos, particularly in this area), hydropower plants, dams, and sensitive infrastructures such as hospitals and schools. The potential infrastructure damage and economic loss from an eruption of the BVF may be significant. A quantitative assessment of exposure to the BVF could help Laotian communities and governments and their partners in planning for future volcanic crises. The national government recently requested such an assessment as part of a geohazards collaboration between the government of Laos and the Earth Observatory of Singapore, through the CCOP (Coordinating Committee for Geoscience Programmes in East and Southeast Asia), in 2019.

Our assessment consists of four parts: i) Identification of vents and lava flows younger than (or indistinguishable in age from the meteorite impact) (~ 790 ka—present); ii) Spatial-density analysis of volcanic vents to inform future vent opening location probability; iii) Probabilistic modelling of lava flow inundation using MOLASSES (Modular Lava Simulation Software for Earth Sciences) from Gallant et al. (2018); and iv) Exposure assessment by combining the flow inundation map with population, infrastructure, and land cover data.

Apart from effusive eruptions, also volcanic explosive activity can occur on the BVF. Here, we qualitatively assess and discuss one type of explosive activity and its hazards, which is phreatomagmatic activity, however, we acknowledge that other types of explosive activity (e.g. Strombolian) might have or may occur on the BVF, but their impacts are generally constrained to areas closer to the vent than phreatomagmatic explosions. Larger explosive eruptions (e.g. sub-Plinian, Plinian) are unlikely to have occurred here. Phreatomagmatic eruptions generally occur when rising magma interact with external water (surface or groundwater) or when moving lava encounters water-saturated sediments. There are still research gaps about this eruption style, but phreatomagmatic activity is thought to be complex, long-lived and capable of producing a diverse range of near-vent hazards (Németh and Kósik, 2020). For example, maar-forming eruptions can produce pyroclastic surges that can travel for kilometers and pass topographic obstacles (e.g. Giordano et al., 2002), and since the source of these eruptions can migrate laterally over time, the potential impacted area can also increase (Graettinger and Bearden 2021). While we have no direct evidence in the products (potentially because of a lack of field data), we consider explosive magma-water interaction as a possible hazardous phenomena on the BVF, given the large amount of water available (surface and groundwater) and the relatively large areal extend from where the magma can rise to or close to the surface.

Finally, we estimate a maximum BVF Average Recurrence Interval (ARI) for eruptions, based on the geological and geochronological information available; because of a lack of data, it is likely that the ARI is shorter.

Geographic and socio-economic aspects of the Bolaven plateau and surroundings

The Bolaven plateau extends over about 6000 km² in the southern portion of the Mekong river basin, in southern Laos (Fig. 1). Champasack province covers the majority of the plateau and shares its northern border with Salavan and Xekong provinces and its southern border with Attapeu province. Most of the plateau has low relief, and summit elevations range between approximately 1000 and 1350 masl. The unique combination of gentle topography, temperature, rainfall, solar radiation, and fertility of the basaltic soils makes the Bolaven plateau the most ideal place for coffee production in Laos (Tulet, 2007; Toro 2012). Rubber, tea, cardamom and corn are also cultivated in the region, but to a lesser degree (Delang et al., 2013). Agriculture is the primary source of income for residents of the plateau and represents about 15% of the nation’s GDP (Delang et al. 2013; Applegate 2016; World Bank 2019).

Hydroelectric power generation is another important sector of Laos’ economy. It accounts for ~ 85% of power production and supply in Laos (Department of Energy Policy and Planning and Ministry of Energy and Mines 2020). The Bolaven plateau has particularly favourable conditions for hydropower production, for the presence of large rivers and elevation gain. About 10% of Laos’ operational dams are on the plateau (Open Development Laos 2016). Of relevance to the aims of our study is the sudden failure of one dam on the southern portion of the plateau, which occurred in 2018, causing floods and fatalities (Latrubesse et al. 2020). This study from Latrubesse and co-authors highlights the general difficulties in this region in assessing potential hazards, mostly due to lack of resources and adoption of proper development strategies (e.g. use of inappropriate material to build dams).

Pakse, the largest city and capital of Champasak province lies on the east bank of the Mekong River, about 15 km from the western edge of the plateau. Smaller cities in the area are Paksong, which sits atop the plateau near the summit of the BVF, and Salavan, about 30 km north of the plateau, on the Xedon River, a tributary of the Mekong (Fig. 2). Numerous villages dot the volcanic landscape. All told, the total population living on and around the plateau is nearly 700,000 (WorldPop 2020). Pakse is an important commercial and transport hub connecting Laos with neighbouring Vietnam, Thailand and Cambodia, and is also a major tourism centre. Paksong and Salavan are key assets for the production and export of coffee, but are also growing centers for tourists interested in exploration of historical and natural sites. Both tourism and agriculture can be impacted by natural hazards. According to a recent study (Japan International Cooperation Agency 2015), Laos experienced, between 1983 and 2012, floods, cyclones, landslides and earthquakes, but only floods and cyclones caused significant damage. In this study, the Japanese agency defines a hazard as “significant” when it causes damage > 1% of GDP and deaths > 100. Surprisingly, the report states that there are no volcanoes in Laos. A general lack of awareness of volcanic hazard may, in fact, explain the lack of national volcanic hazard monitoring and data. While there are several national agencies (e.g. National Disaster Management Office, Water Resources and Environment Agency, Department of Meteorology and Hydrology) in charge of response to floods, cyclones, landslides and earthquakes, there is some ambiguity concerning response to volcanic eruptions. We are keen to help remedy this situation by providing quantitative information that could help communities and local and national authorities improve their ability to respond to future volcanic activity.

Fig. 2

Simplified geological map of the BVF. Volcanic units younger (orange) than the inferred Australasian impact crater cover the majority of the field, likely burying most older units (purple where exposed). Some of the youngest flows are indicated on the map. All mapped vents (black polygons) are scoria cones. Main cities are outlined in white. The Mekong river is flowing southward. Basemap is the shaded relief obtained from the SRTM. Yellow stars indicate locations shown in Fig. 3

Geology of the Bolaven Volcanic Field

The BVF consists of basaltic lavas that were erupted throughout the past ~ 16 Ma (Sieh et al. 2019), through thick flat-lying Mesozoic non-marine clastic sedimentary rocks. Its summit lies on the western portion of the plateau, and flows have spilled over its northern, western, eastern and southern flanks. The flows have spread over about 5000 km² (Fig. 2). Most of the flows are tholeiitic in composition, but many are alkali basalt (Sanematsu et al. 2011; Sieh et al. 2019), and can be rarely directly linked to their sources. Many flows are < 100 ka old, and the youngest flow complex (southern flow) may be barely prehistoric (< 10 ka), judging by the immaturity of virgin forests that cover it (Fig. 3d). The southern flow complex, which erupted from a source on the plateau about 13 km south of Paksong, includes flows that extend southward about 45 km, over the flank of the plateau and onto the adjacent plain (Fig. 2). The bulk volume of the BVF is ~ 900 km³ (Sieh et al. 2019), which encompasses 76 scoria cones. Overall lava flow thicknesses on the BVF range from about 350 m at its summit area to just a few meters on its perimeters.

Fig. 3

Representative examples of lava flows and scoria cones found on the BVF. Approximate ages are provided, based on absolute dating at nearby locations. a 34 ka old lava flow; b Basalt boulders 215 ± 41 ka old, surrounded by saprolitised basalt; c Contact between a flow younger than 40 ka (top) and an older flow, ~ 150 ka (bottom); d Satellite image of a relatively young lava flow (< 40 ka) next to an older one (~ 200 ka), separated by a river (likely displaced by the younger flow). Also note the different colour in vegetation, with the greener vegetation (more mature) for the older flow; e Satellite image (DEM in the inset) of two uneroded, young scoria cones of the young North Flow Complex. Exact locations are shown in Fig. 2 (yellow stars)

Other studies of the Bolaven plateau have focused on mineral resources, such as bauxite, gold and copper (e.g. Sanematsu et al. 2011; Phommakaysone 2012; Long et al. 2019). The Japan International Cooperation Agency (2008) created a geological map of the Bolaven plateau to assist the Lao PDR Ministry of Energy and Mines in identifying zones with high potential for mineral extraction. The JICA report includes geochemical and petrological data from the BVF lavas (pyroxene basalts, andesites, olivine basalts, and alkali basalts) and twelve K–Ar/Ar–Ar dates (Neogene to late Quaternary). Sanematsu et al. (2011) investigated laterite formation on the BVF basalts. Their geochemical and chronological work led them to divide the basalts into three groups (small-volume alkali basalts, 15.7 Ma; large-volume olivine tholeiites, 1.2 Ma; and quartz and olivine tholeiites, younger than 0.5 ± 0.2 Ma). Their five ⁴⁰Ar-³⁹Ar dates are in agreement with Sieh et al. (2019), who provide a far more extensive dataset for the BVF that sheds much more light on its complex volcanic history. In particular, the young ages found for some of the flows (< 40 ka), besides the amount of data recently made available, motivated us to investigate the BVF further, with a focus on the more recent volcanic history for the first time.

Methods

Mapping techniques

We present an updated version of the BVF geological map published by Sieh et al. (2019) (simplified map in Fig. 2 and detailed map in a data repository at https://researchdata.ntu.edu.sg/dataset.xhtml?persistentId=doi:10.21979/N9/HQDXRQ&faces-redirect=true). The updated map includes new field observations, 30 new ⁴⁰Ar-³⁹Ar dates of lava flows, and new interpretations based on digital imagery from ©Google Earth Pro and a geospatial analysis using Esri® ArcMap 10.7.1 (Fig. 3). Field observations included: 1) The use of vegetation as a proxy to determine the relative age of the very youngest flows (i.e. an immature virgin forest implies a very young age), before performing ⁴⁰Ar-³⁹Ar dating; 2) Geomorphological evidence (e.g. differential erosion of neighbouring flows indicating different relative age, or evidence of rivers displaced by lava flows); and 3) Different macro-scale flow textures (e.g. colour, relative abundance of minerals and mineral type, level of weathering). Satellite imagery was mostly used in combination with field observations. ArcMap 10.7.1 was used to create topographic contours from a Digital Elevation Model (Shuttle Radar Topography Mission – SRTM, 30-m resolution), in order to identify, where possible, different flows (even within the same flow complex) and topographic features associated with them (e.g. scoria cones). Scoria cones were mapped based on (i) height (at least a few tens of meters high), (ii) morphology, which also provided clues on their relative age (i.e. conical, horseshoe-like, flat-topped, with the first and the last being relatively younger and older respectively), and (iii) location respect to the nearby 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, Figure S2). Besides scoria cones, other positive-relief topographical features were present, and were mapped as “mounds”. A common characteristic of these mounds is the relatively low H/W ratio (height less than 10 m, occasionally up to 30 m; width from a few hundred meters up to a few kilometers). We interpret them as lava flow structures (e.g. tumuli, lava rises, hornitos, squeeze-ups; e.g. Nemeth et al. 2003; Murcia et al. 2014; Mishra et al. 2019), but we do not exclude the possibility that some may represent volcanic vents, and/or artefacts created by tall, dense vegetation during the generation of the DEM. Also some volcanic fissures were inferred from the SRTM, based on topography and geology.

Kernel density estimation: definition, use, and applications in volcanology

A Kernel Density Estimation (KDE) is a nonparametric statistical method used to represent the distribution of physical features (points or lines) through space. Over the last decade it has been widely used in geoscience (e.g. Ramanna and Dodagoudar 2012; Colkesen et al. 2016). For volcanic hazard assessments, the KDE method is used to predict the location of future vents, based on the location of vents formed in the past (Connor et al. 2012; Bartolini et al. 2013; Gallant et al. 2018). This process assumes that new vents will form through the same magmatic/tectonic mechanisms that formed the old vents (Connor et al. 2019).

Elements of a KDE are the kernel function and the bandwidth. The kernel function defines the probability of future vent formation at locations within a certain region, and can be of different statistical types (e.g. Gaussian, Cauchy, Epanechnikov, Triangular, Uniform, Triweight, Quartic); the bandwidth is the search radius within which the density is calculated. The latter, in particular, is more sensitive to the resulting output relative to the type of kernel function (Martin et al. 2004), and therefore requires careful selection. We performed the KDE with ArcMap 10.7.1, which uses a Quartic kernel function (Silverman 1998) and provides a default bandwidth (defined as search radius in the software) of ~ 7.4 km, computed specifically for the input dataset (distribution and number of scoria cones across the BVF) using a spatial variant of Silverman's Rule of Thumb (Silverman 1998). This method, compared to more classic spatial density estimations, is weighted for spatial outliers, which pertains to some vents on the BVF (Fig. 4).

Fig. 4

Vent spatial-density map of the Bolaven Volcanic Field. The KDE utilized only scoria cones younger than the inferred meteorite impact (< 790 ka)

As input data for the KDE on the BVF, we selected vents based upon two main criteria: 1) Our confidence that the source is indeed a vent; and 2) The vent age is younger or equal to 790 ka (the age of the Australasian impact crater). For the first criterion, we binned eruptive vents on the BVF into three types: Scoria cones, mounds and fissures. Only scoria cones (n = 76) were considered for the estimation of the vent spatial distribution on the BVF. We ignored the many low-relief mounds and fissures, because we know from ground-truthing some of these that the dense and tall vegetation common in the region is the cause of some of these features in the DEM. Moreover, the relative vertical accuracy of the SRTM is around 6-m at these latitudes (Brown et al. 2005), about the same as the height of these features. We constrained vents to those younger than the impact crater age as a reference for the volcanic hazard assessment of the BVF for two main reasons: i) The Australasian meteorite impact represented an important event in the volcanic history of the Bolaven plateau. Although volcanism began long before the impact, there seems to be a spatial relationship between the impact site and the subsequent scoria cones (Fig. 4). We do not rule out the hypothesis that the event may have affected the rate of melt production (discussed later), as has been inferred for other large impacts (Jones 2005). ii) Absolute ages were not available for all lava flows, and this, together with the limitation in mapping remote tropical regions (dense vegetation, high erosion rate, little product exposure) limited our understanding of stratigraphic relations between different flows, and their sources, with a similar state of preservation but different age (e.g. 100 ka vs 200 ka), while the 790 ka meteorite marker better constrained relative ages of lava flows. iii) The chosen ~ 790 ka time window is in line with previous volcanic hazard assessments for similar volcanic fields (e.g. Connor et al. 2012; Gallant et al. 2018). Details on the assignment of scoria cone ages appear in 17 Section.

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

Lava flow simulations

MOLASSES

We use the lava flow inundation model MOLASSES (Gallant et al. 2018) because it requires relatively few input data (lava flow thickness and volume), compared to more complex models (e.g. DOWNFLOW, MAGFLOW, Q-LAVHA) where less easily available parameters such as effusion rate, eruption duration, lava flow temperature and composition are also required (e.g. Favalli 2005; Cappello et al. 2011; Mossoux et al. 2016; Richter et al. 2016). Although MOLASSES does not replicate all the complexities of lava flows in terms of their mechanisms of emplacement, which depend on their physical and chemical properties, it has been proven to effectively replicate flow geometries based on the topography (Gallant et al. 2018; Tsang et al. 2020). Also, MOLASSES has been successfully validated through experimental benchmarking (Dietterich et al. 2017) and comparison with lavas flow emplacement at Tolbachik volcano, Kamchatka (Kubanek et al. 2015). In addition, this software (and earlier versions) have been tested for probabilistic-based approaches in volcanic hazard assessments (Connor et al. 2012; Gallant et al. 2018), which was the intent for our study on the BVF. The basic principle of MOLASSES is that an initial limited volume is erupted, with the volume then distributed to adjacent cells, based on distance from the source cell and/or difference in elevation; this process continues until the volume is exhausted (Connor et al. 2012).

Defining eruption source parameters

The range of eruption source parameters used for our lava flow simulations are listed in Table 1, and include thickness, volume (total volume to be erupted) and pulse volume (maximum volume for each pulse). In order to choose adequate parameters, we used a combination of field data and remote sensing data from the BVF (where measurable) and available data from analogue volcanic fields. Data from the BVF include lava flow thickness, length, surface area and volume. Thickness was measured at a few locations in the field (lava flow edges or front) and through remote sensing methods (through construction of structure contours on ArcMap) at the same lava flows measured in the field, to obtain maximum, minimum and mean values. Lengths and surface areas of lava flows were measured through ArcMap, and volumes of these flows were calculated by multiplying surface area and thickness. We lack a complete record of these parameters for the BVF, because the dense vegetation impedes access, our finite amount of time in the field, and the rarity of exposure of thick flow stacks and their bases. We chose the following volcanic fields as analogues for the BVF, based on their similar geological setting (intraplate distributed volcanism), their basaltic compositions, and the availability of data on lava-flow length, thickness and volume: The Shamiram plateau, Armenia (Connor et al. 2012), the eastern Snake River Plain, USA (Gallant et al. 2018), the Northern Harrat Rahat volcanic field, Saudi Arabia (Murcia et al. 2014), and the Auckland volcanic field, New Zealand (Kereszturi et al. 2012).

Table 1 Eruption source parameters used for MOLASSES. These are averaged values across values obtained from direct (i.e. field survey) and indirect (i.e. remote sensing) measurements on the BVF, and analogue volcanic fields (see Defining eruption source parameters Section for details)

Model inputs

Other inputs for MOLASSES included a Digital Elevation Model (DEM) of the area and previous vent locations. We used the Shuttle Radar Topography Mission (SRTM) DEM, which is available at 1 arc-second resolution (30 m × 30 m at the BVF), and which we cropped to a total extent of ~ 12,000 km² (102 km × 118 km grid) to cover the BVF. The locations of the new vents were stochastically sampled according to the vent spatial density map (Fig. 4). One lava flow was simulated for each vent, with the eruption volume and average lava flow thickness stochastically sampled from within the range considered (Table 1). No correlation was assumed between vent location and eruption volume or lava flow thickness. Volume and thickness were sampled according to a log-normal law, in order to preferentially sample values in the lower/medium range and occasionally in the higher range. This allowed the software to provide lava flow outputs in line with the observations at the BVF and other similar volcanic fields (Connor et al. 2012; Kereszturi et al. 2012; Murcia et al. 2014; Gallant et al. 2018).

We chose to simulate 10,000 individual lava flows, in accordance with other similar probabilistic volcanic hazard assessments (Connor et al. 2012; Gallant et al. 2018). We first ran a sensitivity test to compare the 30-m resolution SRTM with a resampled one at 90-m resolution (100 simulations each), to evaluate the relationship between output resolution and computation time. The computing time was lower for the 90-m resolution DEM (~ 3 h vs ~ 5 h), however, we noticed a larger extent of the lava flows for the 90-m DEM (overall inundated area ~ 20% larger). This is because the 90-m SRTM does not capture minor topographic features capable of stopping or deflecting simulated lava flows, particularly in relatively flat areas. Previous studies using MOLASSES, or former versions of MOLASSES, with a probabilistic-based approach (Connor et al. 2012; Gallant et al. 2018), utilised a 90-m DEM. Connor et al. (2012) in particular highlighted the necessity of using higher resolutions to capture minor topographic obstacles and achieve more reliable results, although the same authors considered the 90-m DEM adequate enough for their study. In areas with relatively flat topography and a paucity of prominent topographic features such as ravines, canyons, mounds and channels, the relevance of the Digital Elevation Model used would be minor. Some studies have used better resolutions (i.e. as high as 1-m/pixel) to conduct volcanic hazard assessments (Capra et al. 2011; Becerril et al. 2014; Deng et al. 2019; Tsang et al. 2020), however, these studies were mostly deterministic and/or conducted on individual volcanic edifices, rather than large volcanic fields, hence involving a lighter computation. Using resolutions finer than 30-m for the BVF may provide better results, however, two main issues would be the availability of such DEMs and the computation time required to run 10,000 simulations. Therefore, based on these considerations we decided to proceed with the 30-m resolution DEM.

Model outputs

MOLASSES outputs a hit intensity map, on which the intensity is the number of times that each grid cell is inundated (hit) by lava flows. By dividing the number of hits in each grid cell by the total number of simulations performed, we obtained a conditional probability of inundation for each grid cell. A conditional probability assumes that an eruption has occurred somewhere in the BVF.

We defined three hazard zones, coloured yellow, orange and red, based on the probability of inundation, at the 90^th (1–78 hits), 50^th (78–390 hits) and 10^th (390–779 hits) percentiles. This choice was somewhat arbitrary but taken to reflect the distributions of hits across the BVF, into areas with relatively low, medium and high probability of inundation.

Population, infrastructure and land cover exposure

We considered population, power lines, power stations, dams, roads, and land cover (i.e. forested and vegetated areas, croplands, built-up areas and water bodies) as critical elements to consider for the BVF because of the intense use of land for agriculture and hydropower production and supply. Other aspects of land use, such as building type and purpose, were not considered here due to the lack and/or reliability of data, or because the exposure to future lava flows was considered very low (e.g. Pakse international Airport).

Population

For most population datasets, the population count (or distribution) depends on a series of factors, such as the availability of census information, the disaggregation methods used, and the spatial resolution. In many cases the available results are estimates only (e.g. Freire et al. 2016; Wardrop et al. 2018; Zhang et al. 2018). The BVF is located in the least economically developed country in Asia (Pink 2016) and one of the least economically developed countries in the world (Delang and Toro, 2011). Moreover, the region of the BVF is predominantly rural; these characteristics make population estimates less reliable when geospatial elements such as nightlights are used as an indication of population density and size (Small et al. 2005). Therefore, we decided to refer to three different, widely used, free and relatively up-to-date population datasets, with different spatial resolutions, in order to provide an indication of the uncertainty in our estimate of people exposed to a future eruption. We used WorldPop2020 (Bondarenko et al. 2020), Global Human Settlement Population Grid 2019 (GHS-POP, Schiavina et al. 2019) and LandScan2019™ (Rose et al. 2019). All these datasets present advantages and disadvantages. WorldPop integrates census data (where available) and other geospatial datasets (e.g. settlement location and extent, land cover, roads, building maps, satellite nightlights, vegetation, topography) to predict population density at a spatial resolution of ~ 100 m. Its advantages are the relatively high resolution and data availability for individual countries. Its disadvantages include the lack of reliable data in countries that have not had a census for a long time (last census in Laos was in 2015) and limited reliability of data in rural areas. GHS-POP also combines census data and settlement information (from Landsat satellite data for target periods), to produce population estimates at a resolution of 250 m; it has a lower resolution than WorldPop, and the last version was released in 2019. However, it has been proven to provide reliable results when compared with other population datasets of lower resolution (Freire et al. 2016). LandScan uses an automated model to integrate sub-national level census counts for each country and geospatial datasets, including land cover, roads, slopes, urban areas, village locations, and high-resolution imagery analysis; the resultant population count is based on a 24-h average. Advantages of this dataset include these: (i) Weighting of population count for each country, based on socio-economic and cultural understanding of the area; and (ii) release of an updated version every year. Some of its limits are its relatively coarse resolution (~ 1 km) and the fact that data seem to be less reliable at the transition between urban and rural areas (Calka and Bielecka 2019).

Infrastructure

We used datasets from different sources for different infrastructure types, based on data accessibility, availability, and date of release. We used Open Street Map (OpenStreetMap Foundation (OSMF) & Contributors 2016) for roads; Open Development Laos (2016) for dams; and World Bank (2019) for power lines and power stations. Following the Open Street Map classification, roads were further subdivided into Type-1 and Type-2 (Table 2). The former consists of roads for public use (e.g. country’s roads, links between cities), the latter includes roads for private or semi-private use (e.g. access to housing, industrial or agricultural use). Note that for the purpose of this paper, roads for exclusive pedestrian use were not included.

Table 2 Road types considered for this work, based on OpenStreetMap definitions

Land cover

For land cover, we used Copernicus Global Land Service (CGLS-LC100 2019), which has a spatial resolution of ~ 100 m. We calculated the total area for each land class based on the cover fraction classification, where each pixel represents a fraction (0–100%) for a particular land class. The classes considered were forest, cropland, built-up and water (seasonal and permanent); other types of natural vegetation were included in the forest class.

All these datasets were overlaid on the probabilistic inundation map and processed through ArcMap 10.7.1, in order to quantify the exposure in terms of numbers of people, count of dams and power stations, length of roads and power lines, and surface area of land cover for each class, within each defined lava flow hazard zone.

Results

Field data

For the purpose of this work, here we concentrate only on post-impact volcanic products and their sources (790 ka or younger). More details appear on the interactive detailed map in the Data Repository (https://researchdata.ntu.edu.sg/dataset.xhtml?persistentId=doi:10.21979/N9/HQDXRQ&faces-redirect=true).

Lava flows

Post-meteorite impact lava flows vary in length between a few hundred meters and a few tens of kilometers (up to ~ 50 km, northern flow complex), for a total areal extent of about 3900 km². Thickness of individual flows was estimated based on observations from locations in the field (see examples in Fig. 3) and additional geospatial analysis through ArcMap, to obtain a range of average thickness between about 6 and 14 m. Volumes were calculated based on the surface area of the mapped flows and thickness information, to obtain estimates in the range ~ 0.06—3.4 km³ (6 × 10⁷- 3.4 × 10⁹ m³). Note that thickness and volume were calculated for eight relatively young flows (≤ 120 ka), chosen to be representative of the latest activity of the field, and because they presented advantages in terms of field exposure, preservation and accessibility, compared to other post-impact lava flows, as highlighted by Sieh et al. (2019), with intense laterization and saprolitization for most basalts on the field (Fig. 3).

Scoria cones

Among the 67 available dates from the BVF, 47 are from post-meteorite impact lavas (Fig. 2 and Table S1, supplementary material). Two of these are from lava within scoria cones. Geochronological and geochemical data from scoria cones are sparse, because they are commonly in thick, untracked jungle and their clastic deposits are more highly weathered than lava flows (only one scoria cone was successfully dated, ~ 1.5 km E of Paksong, 200 ± 9 ka). Therefore, most scoria cones were either assigned an absolute age, based on ⁴⁰Ar-³⁹Ar dates of lava flows that were inferred to have erupted from that cone (e.g. Fig. 3e), or assigned an age interval (if an absolute age was not available for any of the flows that emanated from beneath that cone). Out of the 97 scoria cones mapped (Fig. 2), 76 were inferred to be younger than 790 ka, whereas 21 were either older than 790 ka or indeterminate in age/age-interval.

Here we assume that each vent was produced by a single eruptive event (following the definition of vents vs events as in Gallant et al. (2018). This assumption was based on the lack of field information for most scoria cones (lack of fresh scoria cones deposits and/or lack of accessibility), resulting in missing information usually needed to link an eruptive source to its deposits; for example, use of geochemical compositions and absolute ages of scoria cones, to match (or not match) them with nearby lava flows and other scoria cones (more discussion around the issue of vents vs events is reported later in the text). Understanding the direct stratigraphic relationships between scoria cones and their products (i.e. lava flows) was rarely possible. Satellite images could be used to map lava flows from source to end with confidence in only two instances in the BVF (the southern and southwestern lava flow complexes: Fig. 2 and Figure S2). Additional remote sensing methods (e.g. use of topographic contours to understand possible flow direction and provenance) were needed to link lava flow to all of the other scoria cones, but we acknowledge that this comes with uncertainties (i.e. vegetation affecting the accuracy of topographic contours, creating artefacts in the SRTM).

Vent spatial density and probabilistic inundation map of the BVF

The approximately 30% of scoria cones within the highest density fields (5^th to 33^th percentile) form a N-S band across the NW portion of the plateau (Fig. 4). Some isolated clusters of lower density are to the NE (~ 15 km south of Salavan), to the W (~ 10 km north-east of Pakse), and about 35 km N of Pakse. Curiously, the area with the highest concentration of vents is nearly coincident with the inferred location of the Australasian impact crater (Fig. 4).

As one might expect, the highest probabilities of inundation by future eruptions of lava are in the region of greatest density of existing vents. The area with the highest probability of inundation (red zone: 3.89–7.79%) is nearly coincident with the area with the highest density of vents, whereas the lowest-probability area (yellow zone: 0.01–0.78%) has the lowest vent density. The lowest conditional probability of lava flow inundation (yellow zone) has the largest total surface area, at 3508 km², while the orange zone and red zone cover 1804 km² and 285 km², respectively. A total area of approximately 5600 km² is at threat from potential inundation by lava flows in a future eruption. For any individual eruption, the area affected ranges from ~ 2 km² to ~ 740 km².

Exposure

Population

The three population datasets yield a range of ~ 274,000 to ~ 358,000 people within the potentially inundated area. Most (~ 189,000 to ~ 219,000) live in the yellow zone. About 44,000 to ~ 117,000 live in the orange zone and ~ 12,000 to 26,000 live in the red zone. All three of the large cities are within the lowest of the three zones. Stated in percentages and using the WorldPop2020 and GHS-POP 2019 databases, 60% live in the yellow zone, 33% in the orange zone, and 7% in the red zone. The values using LandScan2019, are 80% in the yellow zone, 16% in the orange zone, and 4% in the red zone (Table 3).

Table 3 Population exposure results from different datasets

Infrastructure

There are 2678 km of roads within the three potentially inundated areas (Fig. 5). Of these, 1479 km (55%) are within the yellow zone, 950 km (35%) within the orange zone, and 250 km (10%) are within the red zone. In terms of road type, the balance is slightly in favour of Type-1 roads for yellow and orange zones (55 and 59% respectively), and roughly the same for the red zone (Type-1 = 49%, Type-2 = 51%). This reflects the underlying proportion of roads, which are approximately 56% Type-1 and 44% Type-2 across our study area (see also Fig. 5).

Fig. 5

Probabilistic inundation map of the BVF and exposed key infrastructure for the three hazard zones. Percentage of power lines within each zone as a pie chart in upper right, and percentage of road types for each zone as pie charts at the bottom

Across the entire hazard area are 416 km of power lines. Of these, 259 km (62%) are in the yellow zone, 120 km in the orange zone (29%), and 37 km (9%) in the red zone. Along these power lines, there are two power stations (both in the yellow zone) and one dam (orange zone). Another dam is located in the orange zone (Table 4, Lat: 15.35°; Long: 106.31°), but not yet connected to the main array of power lines (the Laos government planned to build additional 578 km of power lines in the region, including 5 km to connect this dam to the main network, according to the World Bank, 2019).

Table 4 Exposure results for Dams and Power stations

Land cover

Across all hazard zones (Fig. 6), 4996 km² (89%) of land is covered in forest and sparse vegetation. Cropland covers 514 km² (9%), built-up areas cover 73 km² (1%), and permanent and seasonal water cover 15 km² (< 1%). In particular, the yellow zone encompasses 3060 km² (87%) of forest and sparse vegetation, 387 km² (11%) of cropland, 48 km² (1%) of built-up areas, and 13 km² (> 1%) of water. The orange zone contains 1686 km² (94%) of forest and sparse vegetation, 98 km² (5%) of cropland, 19 km² (1%) of built-up areas, and 1.4 km² (> 1%) of water. The red zone includes 250 km² (88%) of forest and sparse vegetation, 29 km² (10%) of cropland, 5 km² (2%) of built-up areas, and 0.6 km² (< 1%) of water.

Fig. 6

Land cover exposure map. Shades of the same colour type represent fractions (e.g. light pink = 1% of cropland, dark pink = 100% of cropland, for that particular pixel), that have been considered for the total calculation of surface area (see pie charts) for each land class

Discussion

Volcanic hazard assessments require information on several aspects, such as volcanological setting (e.g. central volcano, volcanic field), eruptive style (e.g. explosive, effusive), hazard type (e.g. ash fallout, pyroclastic flows, lava fountains, lava flows, lahars), and the geographic context of the area assessed (how the hazards can impact the local population and/or their activities). Their robustness often depends strongly upon the availability of data for that particular area (reliable and up-to-date geological, hazard, exposure and vulnerability information). Here we did not conduct any vulnerability assessment due to the lack of information. Below, we discuss the results and limitations of the present research, based on evaluating likely hazard, exposure and potential impacts in the BVF area in case of renewed volcanism, the first such assessment in northern SE Asia for a volcanic field.

Scoria cones distribution and significance

Volcanic fields around the world show a variety of spatial features of their eruptive centres, from clustered monogenetic cones (e.g. Connor et al. 2012; Kereszturi et al. 2012; Gallant et al. 2018) to cone alignments along fissures (e.g. Belousov et al. 2015; Kubanek et al. 2015; Pedersen et al. 2017). Understanding the distribution of eruptive centres within a volcanic field may help understand the likely location of future activity at that field, hence allowing more accurate hazard and exposure assessments. The location of the post-impact scoria cones on the BVF, hence the area with the highest probability of future vent opening, seems related to the impact crater structure rather than to the local tectonic stress.

From a regional tectonic perspective, several works report the current state of stress field for northern Southeast Asia, including Laos (Pailoplee et al. 2009; Pailoplee and Charusiri 2017; Meyers et al. 2018). Particularly, these works report seismogenic faults in central and northern Laos, which can be divided into two groups, SW-NE and SE-NW oriented faults (Pailoplee and Charusiri, 2017), mostly in a strike-slip regime driven by the Sagaing Fault and Red River Fault (Meyers et al. 2018). None of this work shows local faults and their orientation in southern Laos, where the BVF is located. From a more local tectonic perspective, a report from the Ministry of Energy and Mines (Japan International Cooperation Agency 2008) and mapping efforts from this work (https://researchdata.ntu.edu.sg/dataset.xhtml?persistentId=doi:10.21979/N9/HQDXRQ&faces-redirect=true), highlight the presence of several tectonic structures on and around the Bolaven plateau (a syncline crossing the plateau, a syncline/anticline pair just NE of the plateau, and several reverse faults to the E of the plateau), whose σ1 is oriented NE-SW. Field evidence and Ar–Ar dating of a nearby lava flow indicate that these structures (particularly the syncline/anticline pair NE of the plateau) are younger than 2.9 Ma, and they likely represent the latest deformation pattern in southern Laos. If there was a direct link between this local stress field and the pattern of scoria cones, then it would be parallel to the σ1 (Nakamura, 1977), which is not the case here, given the N-S orientation of the scoria cones.

On the other hand, the N-S aligned post-impact scoria cones on the BVF, besides being encompassed within the inferred location of the Australasian impact crater, have an inferred age lower than the impact itself, with pre-impact scoria cones showing a more widespread distribution across the field (Fig. 4). This may suggest an existing link between the meteorite impact and the post-impact volcanism, as hypothesized for other large-scale meteorite impacts (Jones 2005), or in part may reflect a bias from the lack of exposure of older scoria cones (buried or eroded), potentially hiding older (or unknown) regimes.

Although additional data are needed (including a complete dataset of seismic tomography, 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 location of future volcanic activity on the BVF, which largely coincides with the location of the largest known young meteorite impact on Earth.

Exposure analysis

In this section we discuss exposure and impacts on populations, infrastructures and landcover, based on our lava flow simulations and mapping.

The population exposure results show that there is some variation among all datasets. In particular, the highest variation is observed for LandScan, where the total count is between 13 and 24% smaller than that for WorldPop and GHS-POP (Table 3) respectively, and also different for individual hazard zones (20% higher for yellow zone, 17% smaller for orange zone, and 3% smaller for red zone). GHS-POP has the highest total population count (~ 358,000), but if we consider the fractions for each hazard zone, WorldPop and GHS-POP share the same % of people. All three datasets are based on census data and geospatial information. Decadal censuses are conducted in Laos, with the last conducted in 2015, hence relatively recent. A major difference among the considered datasets is the spatial resolution, with LandScan having the coarsest one. In addition, as anticipated in the Methods Section, LandScan has been found to underestimate population counts at the transition between urban and rural areas (Calka and Bielecka 2019) and the lower LandScan population count on the BVF may reflect the presence of sparse villages in the region, surrounded by forests and croplands (Fig. 6). WorldPop and GHS-POP may present some limitations as well, however, the results are relatively consistent when total counts and fractions across the different hazard zones are compared. For the purpose of this study, the application of largely used, and freely available population datasets allowed us to detect the minimum number of people potentially exposed to lava flows in case of an eruption on the BVF.

Generally, casualties are less likely from lava flows than from pyroclastic flows, ash fallout, ballistics, or lahars (Brown et al. 2017). However, lava flows can significantly impact settlements and livelihoods, and rehabilitation/recovery can take years and involve high costs for the local administrations (Williams and Moore 1983; Baxter et al. 2002; Jenkins et al. 2017; Tsang and Lindsay 2020; Meredith et al. 2022). An exhaustive review on impacts from lava flows is provided by Harris (2015). He considered six main categories potentially impacted from lava flows: 1) humans and animals; 2) buildings; 3) communication, utilities and machinery; 4) agriculture, habitat and natural vegetation; 5) economic activity; and 6) social fabric. All these categories would have a different impact based on factors such as the characteristics of the flow (e.g. viscosity, advancement rate, 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 strategies (e.g. redirection of lava flows to less exposed/vulnerable areas). Education of the population to the likely hazards posed by lava flows, and their long-term impacts can help communities adapt and recover (Tsang and Lindsay 2020).

Cities on and around the BVF have calculated conditional probability of inundation that vary from 0.01% to 7.79% across the field (Fig. 5). Pakse is the largest and most populated city in southern Laos, and ⁴⁰Ar-³⁹Ar dating shows that the majority of the city is built upon a lava flow ~ 180 ka old, and a flow likely < 40 ka is nearby. This low spatial density of past flows puts Pakse within the least exposed zone, our Yellow Zone (Fig. 5). It has about a 0.05% likelihood of being reached by lava flows during any future eruption. That is to say, there is one chance in 2,000 that any future flow will enter the city. Paksong, the second largest city in the area, sits within the Red Zone. Although the flow that it sits upon is about the same age as the one underlying Pakse (~ 185 ka), Paksong is within a few kilometers of many young flows that range in age from ~ 75 to ~ 200 ka. The likelihood of a future flow reaching Paksong is ~ 6% (or about one chance in 17), about a hundred times higher than one reaching Pakse. The city of Salavan is north of the BVF and is built on young fluvial sediments rather than lava flows. Our simulations thus show a conditional probability of inundation of 0%. That is not to say that there is absolutely no possibility of inundation by lava, but the likelihood is vanishingly small. Salavan is more likely to be affected by other aspects of an eruption. For example, a flow down the north side of the plateau could reach and dam the Xedon River, flooding parts of the plain around and including the city.

People can be significantly affected by damage to infrastructure. Disruption of power production and distribution could result in loss of access to power for most areas in the Champasak, Salavan, and Xekong provinces. Here we only consider hydropower, for the reasons explained in previous sections, which is produced by the operational dams in the area, and distributed through a network of power lines (Fig. 5). Eruption-related damage to one of the two dams located in the orange zone, which is directly connected to one of the power lines, may create a disruption in the immediate power production and supply. Power lines on and around the BVF cover the entire width of the potentially inundated area, therefore, regardless of where the damage will take place, it may affect other areas of the plateau, potentially also affecting urban sites located outside of this area (e.g. Salavan). In addition, most of the power produced in Laos is exported to Thailand; therefore, international trading can undergo an impact as well. Another key infrastructure that can suffer an impact, affecting people, is roads. This can happen through restricted or blocked access for: 1) delivery of essential goods (e.g. food, medicines) to the areas affected; 2) emergency response or maintenance vehicles, for example to 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 villages who need assistance; 4) trading with nearby countries (Thailand, Cambodia and Vietnam).

Another key asset to consider in exposure assessments is land cover. Our analysis of the potentially inundated area shows an obvious predominance of forest and other types of natural vegetation, followed by croplands, built-up areas and water bodies (Fig. 6). An eruption on the Bolaven plateau affecting them can result in impacts for the local and national economy. Forests, for example, may be ignited by lava flows (e.g. Ainsworth and Kauffman 2009; Harris 2015), particularly in the drier season, and fires can propagate over a large area, potentially damaging infrastructures and natural historical sites. Cropland and built-up land classes, which seem to be linked across the whole field (built-up areas being surrounded or adjacent to areas with high cropland fraction), double from orange zone to red zone (5% to 10% for cropland, and 1% to 2% for built-up areas), despite the red hazard zone being ~ 6 times smaller than the orange hazard zone; this suggests that the high area of the plateau (also the red zone: Fig. 6) has the most favourable conditions for agricultural use. Although cropland was not subdivided into different agricultural classes here (lack of reliable and/or up-to-date data), coffee production on the Bolaven plateau represents about 95% of the total amount produced in Laos (Toro 2012). A future eruption on the Bolaven plateau could therefore affect the socio-economic wellbeing of those in the region, who rely on coffee as their source of income. This in turn can impact the economy of the whole country (Toro 2012), either directly, if coffee plantations are inundated by lava, or indirectly, if roads linking cities or countries are inundated by lava and inaccessible for months/years affecting transport to major cities for coffee processing and export. Official GDP data from Laos are available for the agriculture sector as a whole, but not available for individual sub-sectors such as coffee.

In order to further evaluate the knock-on consequences from an effusive eruption on the BVF, up-to-date information about infrastructure and land cover is needed.

The hypothetical role of external water in explosive volcanic activity on the BVF

Explosive interaction between magma and water (here broadly referred to as phreatomagmatic activity), is a potentially hazardous volcanic phenomenon known to occur in a wide range of environmental settings (e.g. Thorarinsson 1967; Fagents and Thordarson 2007; Wohletz et al. 2013; Verolino et al. 2018, 2019; Dürig et al. 2020a, 2020b), including areas dominated by lava flow effusions (e.g. Lorenz and Haneke 2004; Hamilton et al. 2010; Fitch and Fagents 2020). Although we found no direct evidence of past explosive magma-water interaction on the BVF, the abundance of water in this area may lead to such activity in case of eruption. Below we address three different potential mechanisms leading to such behaviour on the BVF: (i) new vents opening underneath standing water (reservoirs or rivers) (Fig. 7); (ii) rising magma interacting with groundwater; and (iii) lava flows interacting with surface water or water-saturated sediments. The hazards associated with all these mechanisms are similar, including formation of pyroclastic surges, delivery of ash, lapilli and bombs into the atmosphere and to nearby areas. In addition, if an eruption initiates in a large water reservoir, the consequence can be generation of a tsunami (e.g. Waters and Fisher 1971; Vaughan and Webley 2010; Sandri et al. 2012; Nomikou et al. 2014; Deligne et al. 2017). All these hazards can be more dangerous to people than lava flows, but also can cause damage to infrastructure and land cover (Jenkins et al. 2015).

Fig. 7

Probabilistic inundation map with the location of the simulated vents (n = 10,000) and simulated vents that opened within water bodies (n = 31)

Among the 10,000 simulated vents, 0.31% of them (n = 31) opened beneath standing water (Fig. 7). Explosive interactions between standing water and rising magma, although more common in oceans, do occur in lakes (e.g. Németh et al. 2006; Verolino et al. 2018, 2019) and rivers (e.g. Muller and Veyl 1956; Hamilton and Myers 1963; Hackett and Morgan 1988). Should an eruption start on the Bolaven plateau within water, it would likely be phreatomagmatic in style. Paksong is the city with the highest chance of being inundated by lava flows, because of the high concentration of vents nearby; among the submerged simulated vents, two are located less than 500 m from its administrative boundaries, and eight are located ~ 4 km away.

Interaction between rising magma and groundwater on the BVF is also possible. The Bolaven plateau represents one of the rainiest regions in Laos, with a long-term precipitation average of ~ 2500 mm/year (Laos average ~ 1500 mm/year) (Viossanges et al. 2017). This, together with the geological and geomorphological characteristics of the area, make it an exceptional place for groundwater storage, productivity and recharge (Viossanges et al. 2017). Having such availability of groundwater beneath a volcanic field may result in phreatomagmatic activity, in case of renewed volcanism in the area. Explosive interactions between rising magma and groundwater are more likely to take place if they initiate at depths < 1 km below the water table (Valentine and White 2012), resulting in the formation of a maar-diatreme volcano. Borehole data suggest that groundwater of the Bolaven plateau can be found at a depth of a few tens of meters (Viossanges et al. 2017), hence retaining the potential to trigger phreatomagmatic eruptions with centers proximate to the surface, resulting in more hazardous scenarios. The extent of the Bolaven plateau (~ 6000 km²), and its abundance of groundwater, suggests that the probability of this type of volcanic activity is clearly higher than phreatomagmatic activity produced by magma and standing water.

Ultimately, fresh lava flows on the Bolaven plateau can interact with surface water (water reservoirs or rivers) and in some case can result in secondary phreatomagmatic activity (i.e. there is no involvement of new magma, and the explosions can take place even kilometers away from the vent), forming rootless cones (e.g. Boreham et al. 2020; Fitch and Fagents 2020). Secondary phreatomagmatic explosions are more likely if there is interaction between fresh low-viscous lava and water-saturated sediments in an unconfined environment (Zimanowski and Wohletz 2000; Schipper et al. 2011). In most cases they are relatively weak, with erupted products confined to maximum a few ten of meters from the secondary vent (e.g. Fagents and Thordarson 2007; Fitch and Fagents 2020), however, Belousov et al. (2011) showed that they can be equivalent to primary phreatomagmatic eruptions in terms of magnitude and potential hazards (formation of a vertical eruptive columns up to 7 km high, and pyroclastic flows travelled up to 2 km from the explosion site), if they occur on steep slopes. A common scenario of these interactions is a lava flow entering a river, lake, or ocean, without causing an explosion. The lava solidifies and creates new land, resulting in deflection (e.g. Pedersen et al. 2017) or blocks the river, with consequent upstream flooding (Boreham et al. 2020). Approximately 60% of the simulated inundation area of the BVF is covered by post-impact lava, with a good part of it turned into saprolite, and the remaining ~ 40% being Mesozoic (mudstone and sandstone) and Holocene fluvial/lacustrine sediments (loose coarse silt to fine sand). This area encompasses a total of about 15 km² of surface water (Fig. 6) with ~ 1400 km of rivers/streams (Fig. 7), and this estimate does not include the water-saturated sediments, that may be found either in the saprolites/laterites, or in the loose fluvial/lacustrine units. Therefore, this analysis shows that there is potential to trigger secondary phreatomagmatic explosions on the BVF, as well as non-explosive interactions, which can result in damming of a river. Both scenarios can be hazardous for the local communities, either directly (e.g. pyroclastic flows, ejecta, lahars), or indirectly (e.g. flooding and power disruption). However, given the high variability of the physical properties of the sediments involved in the area, the complexity of the topography (flat areas versus steep canyons), and the availability of water which depends on seasons, we cannot say with what frequency in time and space they will occur.

As anticipated at the beginning of this section, we found no evidence of explosive magma-water interactions (i.e. tuff rings, maars) on the BVF, after a qualitative morphologic analysis of the scoria cones. One common feature of tuff rings and maars is their deep craters (i.e. crater base being at the original pre-eruptive surface level, or even lower in the case of maars). The three main vent morphologies found in the field include sub-conical shapes with shallow crater base, horseshoe-like cones and flat-topped ones. All of them have been interpreted as scoria cones, shaped by the eruption that formed them, and/or by erosion. There may be several reasons behind the lack of evidence of phreatomagmatic activity. One reason could be that pyroclastic material can be rapidly removed or mobilised during and after the eruption (e.g. Ferrucci et al. 2005; Romagnoli and Jakobsson 2015; Németh et al. 2021; Simurda et al. 2022), and we would expect this process to be even more prominent in terms of long-term weathering in a wet and warm climate (Németh and Cronin 2007). Another reason could lie in the availability of water, which usually drives the eruption style (phreatomagmatic versus magmatic). This can depend either on local conditions at the scale of a single eruption (depletion of water in the aquifer as the eruption progresses) or regional water table fluctuations due to paleo-climatic changes over an entire eruptive cycle (Németh et al. 2001; Kereszturi et al. 2010; Kshirsagar et al. 2016; Agustín-Flores et al. 2021). Heaney (1991), for example, proved through fossil pollen records that the climatic conditions in the Late Pleistocene were much drier than today for this region. Most of the youngest lava flows on the BVF are latest-Pleistocene, although some flows may well be Holocene in age. At this time, the climate was shifting toward drier conditions until reaching the Last Glacial Maximum (~ 18 ka; Heaney 1991). Therefore, it is likely that at the time of the most recent eruptions on the BVF there was less availability of water than nowadays. Landforms resulting from past phreatomagmatic activity in this area may also be buried under younger lava flows, and their pyroclastic products may be hidden by the dense vegetation, which still makes several areas of the plateau inaccessible. We acknowledge, however, that our discussions are hypothetical in nature. For example, the fact that most magmas have to pass through groundwater or surface water on their way to eruption, but only few of these eruptions are phreatomagmatic, makes them rare events, and highlights that “special” conditions are needed as a trigger mechanism. One factor that seems to favour these eruptions is the pulsing nature of the magma-water interaction, for instance, Doubik and Hill (1999) attributed variations in eruptive style (from dry magmatic to phreatomagmatic) at Tolbachik volcano to conduit-wall collapses, which in turn allowed groundwater to flow downward during stages of low conduit-flow pressure, triggering phreatomagmatic explosions. Another favourable condition is the fractured nature of an aquifer, which allowed more efficient phreatomagmatic explosions at Quaternary monogenetic volcanoes in western Mexico (Agustín-Flores et al. 2021). The subsurface conditions at the Bolaven Volcanic Field are not known well enough to make any of these hypotheses; however, future targeted field campaigns may unveil whether or not phreatomagmatic activity occurred in the past in this area, and from there a quantitative phreatomagmatic hazard assessment could be conducted.

On the other hand, non-explosive magma-water interactions largely occurred on the BVF; one example of river diversion is reported in Fig. 3d.

Average Recurrence Interval (ARI) of the BFV

Evaluating the average recurrence interval at volcanic fields is key in volcanic hazard assessments. Previous authors have estimated ARI’s for volcanic fields by simply dividing a considered timeframe by the number of eruptions that occurred within that period. One common challenge in these assessments, however, is to distinguish eruptive vents (i.e. one vent = one eruption) from eruptive events (i.e. one or more vents related to the same eruption). Runge et al. (2014) and Gallant et al. (2018) attempted this distinction at the Northern Harrat Rahat volcanic field and Eastern Snake River Plain respectively, to provide the minimum number of eruptions. In particular, they used statistical methods based on spatial–temporal relationships of the vents to convert them into events and obtain the recurrence intervals at the selected locations. However, the areas investigated have high levels of exposure (e.g. little or no vegetation), which made it easier to understand the stratigraphic relationships between lava flows and/or vents, and have a well-known volcanic history with absolute ages available for all the vents.

Since the BVF sources are neither well exposed nor thoroughly dated, we approached the challenge of estimating eruption recurrence by considering two sets of data. One set comprises all the post-meteorite impact absolute ages for the BVF lava flows, which differed significantly in age and/or in space. We considered any two neighbouring flows with significantly different absolute ages to represent separate eruptive events. This analysis resulted in 40 eruptions in the last 790 ka, equivalent to an ARI of 19,750 years (Fig. S3). If we apply the same exercise to the last 135 ka, the number of eruptive events is 15, resulting in an ARI of 9,000 years. This difference is likely due to the fact that the latest period of activity on the plateau is more represented in the field, and probably more representative of the latest ARI at Bolaven, hence more realistic. Instead, if we consider the whole post-impact period (i.e. 790 ka), we may lose information from older (buried?) flows, resulting in a longer ARI. The second set, instead, includes the scoria cones for the ARI estimate. If we assume that all 76 cones represent separate events over the past 790 ka, we calculate an ARI of 10,395 years. Although it is possible that some of these scoria cones might have formed during the same eruptive event, it is also very likely that some of the eruptions did not produce any cone, and that we missed detecting some eruptions because their cones are buried or eroded beyond recognition. This estimate is in line with the ARI estimate from the lava flow absolute ages for the last 135 ka.

From a field perspective, there are three flow complexes that are covered by immature first-growth forests. The southern flow complex, reported in Fig. 3d (< 40 ka), has the sparsest such cover; the western flow complex (~ 25 km N of Pakse) also underlies a first-growth immature forest, as does most of the southwestern flow (~ 30 km SSE of Paksong). ⁴⁰Ar-³⁹Ar dates thus far on both of the SW and W flow have not been successful, due to lack or scarcity of radiogenic argon. Dating of a small flow near Pakse and another flow in the northern flow complex (~ 20 km SW of Salavan) were also unsuccessful and for the same reason. The implication is that all these flow complexes are likely very young (the ⁴⁰Ar-³⁹Ar dating method becomes less reliable with decreasing age; Hill et al. 1993). From qualitative observations at some locations in the northern flow complex, we also see a very little developed soil compared to nearby flows dated ~ 100 ka. The radiocarbon dating method would offer a better solution for flows younger than 40 ka on the BVF, however, organic matter is unlikely to be preserved in basaltic lava flows, where temperature exceeds 1000 °C, and has not been found in any of the flows, nor in any paleosol between them. Further geological/ecological studies would be valuable to map and date the flows across the BVF in more detail and provide a more precise ARI. For example, Investigation of vegetation re-growth on lava flows would be a valuable tool to estimate their relative age, as demonstrated in several instances (e.g. Li et al. 2015, 2018), however, this method appears to be more efficient for relatively young lava flows (a few tens or hundreds of years old) (Li et al. 2015, 2018; Albert et al. 2020). The southern and southwestern lava flow complexes on the BVF are relatively young but in the order of thousands of years, however, their vegetation differences compared to the surrounding lava flows make them the only flows on the BVF distinguishable from satellite imagery analysis. In such a tropical climate we would expect very rapid vegetation colonisation, but a quantitative estimate is difficult to achieve. Other factors may influence natural vegetation colonisation in different tropical/sub-tropical climates, such as elevation, slope, local variations of temperature and rainfall, and native soil type, among others (Li et al. 2018). A study on a similarly wet but not as warm climate at Rangitoto volcano (Auckland Volcanic Field, New Zealand) (Wotherspoon and Wotherspoon, 2002) highlights how the vegetation colonisation on lava flows can be strongly controlled by the lava type. In their study, the authors report that about 80% of the lava flows have been colonised by native vegetation since the eruption that occurred ~ 600 years ago. Of these, pahoehoe lavas are mostly covered by vegetation, while most aa lavas are only covered up to 60%. The same authors argue that it may take hundreds of years before all these flows are fully clothed in natural forests. This may be the case for some of the young flows on the BVF, where, although not always possible, pahoehoe and aa lavas were both recognised in the field at several sites (see field notes in the data repository: https://researchdata.ntu.edu.sg/dataset.xhtml?persistentId=doi:10.21979/N9/HQDXRQ&faces-redirect=true). Further on-site investigations will be required in this regard.

Our ARI estimate likely retains high levels of uncertainties, what is certain, however, is that the last eruptive activity on the Bolaven Volcanic Field is more recent than previously thought, as also supported by field evidence.

Conclusions

We conducted the first volcanic hazard and exposure assessment of a volcanic field in northern SE Asia, using field and geochronological data, remote sensing and numerical modelling. We find that future eruptions of lava in the BVF may have socio-economic impacts on the local population and their activities. Inundations of cities and villages, roads, power lines, dams, and coffee plantations are possible but have low likelihood relative to more active volcanic arcs in southern SE Asia. The abundance of groundwater in the weathered and permeable lavas of the plateau mean that explosive phreatomagmatic activity is also plausible and may be more hazardous than lava flows. The maximum average recurrence interval of the BVF is ≤ 10,400 years, much shorter than the time since the 34 ka date of the most recent Ar–Ar-dated eruption. This initial study of the BVF provides the tools and emphasizes the necessity of investigating other understudied volcanic fields in the region, in order to improve awareness of and preparedness for future volcanic crises.