We present a new probabilistic lava flow hazard assessment for the U.S. Department of Energy’s Idaho National Laboratory (INL) nuclear facility that (1) explores the way eruptions are defined and modeled, (2) stochastically samples lava flow parameters from observed values for use in MOLASSES, a lava flow simulator, (3) calculates the likelihood of a new vent opening within the boundaries of INL, (4) determines probabilities of lava flow inundation for INL through Monte Carlo simulation, and (5) couples inundation probabilities with recurrence rates to determine the annual likelihood of lava flow inundation for INL. Results show a 30% probability of partial inundation of the INL given an effusive eruption on the eastern Snake River Plain, with an annual inundation probability of 8.4 × 10−5 to 1.8 × 10−4. An annual probability of 6.2 × 10−5 to 1.2 × 10−4 is estimated for the opening of a new eruptive center within INL boundaries.


The intersection of volcanic hazards and sensitive infrastructure, such as nuclear facilities, can be devastating (Menand et al., 2009). Forecasting long-term volcanic hazards (lava flows, tephra fallout, pyroclastic flows, etc.) is an essential step in mitigating risk, as is upgrading existing forecasts as new information (e.g., new unit ages or renewed activity) and modeling technologies become available. Traditional volcanic hazard assessment methods focused on cataloging activity in a region as a proxy for future activity (e.g., Kuntz and Dalrymple, 1979; Hackett et al., 2002). Modern assessment tools combine geologic history with computational methods to improve forecasts (e.g., Tonini et al., 2015). The results of many contemporary forecasts are simulations that present conditional probabilities of activity—probabilities that are dependent on the assumed occurrence of a future event (Connor et al., 2012). Coupling these conditional probabilities with activity rates allows unconditional probabilities to be resolved (i.e., probabilities in terms of time scales). Unconditional hazard probabilities provide a metric that associates time and potential magnitude with activity, which can be useful for evaluating risk (Cappello et al., 2015).

We present an unconditional probabilistic lava flow hazard assessment for the U.S. Department of Energy’s Idaho National Laboratory (INL, Idaho, USA) on the eastern Snake River Plain (ESRP). The INL covers 10% of the ESRP and contains the highest density of nuclear facilities on Earth (INL, 2016). Previous volcanic hazard assessments for the region (e.g., Kuntz and Dalrymple, 1979; Hackett et al., 2002) cataloged previous eruptions and assigned hazard levels based on proximity to young flows. Hackett et al. (2002) reported annual inundation probabilities of 1 × 10−6 to 4 × 10−7 for the Central Facilities Area, located in the southwest corner of INL. Our assessment utilizes prolific geologic mapping and differs from earlier works by incorporating novel models of ESRP volcanism, a new method of clustering vents into eruptive events, probabilistic selection of input parameters, computational lava flow simulations, and analysis of activity recurrence intervals to report unconditional probabilities of future hazards. It is the first to calculate the probability of vent formation within INL boundaries and consider the likelihood that lava flows will cover part of the INL using Monte Carlo simulation.


The ESRP is a 350 km × 100 km depression that subsided in the wake of the Yellowstone hotspot over the past 10 m.y. (McQuarrie and Rodgers, 1998) (Fig. 1). Bimodal basalt-rhyolite volcanism and sedimentation are the prevailing constructional processes, accompanied by subsidence at multiple scales across the ESRP (McQuarrie and Rodgers, 1998; Wetmore et al., 2009). Total basalt thickness exceeds 1.9 km along the northeast-trending axis of the ESRP and tapers to <30 m at the margins (Kuntz et al., 1992; Shervais et al., 2012). Nearly 95% of the ESRP surface volcanics erupted as basaltic shields, cones, and lava flows during the past 730 k.y.; the remaining consist of older basalts, rhyolite domes, and tuff (Kuntz et al., 1994, 2007).

The >500 mapped volcanic vents of the ESRP form a northeast-southwest–trending band along the hot spot track and include several clusters consisting of tens to hundreds of mostly monogenetic vents (Wetmore et al., 2009). Basaltic activity shows no spatio-temporal trend in distribution (Kuntz et al., 1992). Basalt accumulation for a given locality on the ESRP is uniform, although hiatuses of up to 200 k.y. may occur before accumulation resumes at the original rate (Champion et al., 2002). Recent activity consists of <1–6 km3 basaltic fields that are composed of pāhoehoe with minor near-vent ʻaʻā. The broad spatial distribution of recent activity, median repose intervals of 1000–10,000 yr, and relatively homogeneous plain-wide olivine tholeiite basalt composition suggest that magma generation is rapid and episodic beneath the ESRP, with short residence time in upper mantle and/or crustal reservoirs (Kuntz et al.,1992). Fractionation of olivine tholeiite is responsible for more alkalic compositions on the ESRP (e.g., Craters of the Moon National Monument and Preserve [COM]) (McCurry et al., 2008). Although COM volumetrically dominates Holocene volcanism, it should not heavily contribute to a defining long-term eruption model for the ESRP because the overwhelming majority of activity is monogenetic and compositionally primitive (Kuntz et al., 1992; McCurry et al., 2008).


Basaltic eruptions are complex multi-phase processes that can persist for years to decades with pauses and shifts in the location of activity over time (Cashman and Mangan, 2014). These variations raise important questions about how an eruption is defined, for example: does the continuous activity of the 1983–2018 of Puʻu ʻŌʻō (Hawaiʻi, USA) count as many individual eruptions or as a single eruption (Orr et al., 2015)? Can a single eruption occur simultaneously from multiple vents, such as the 2012–2013 Tobalchik flows (Kamchatka, Russia), or does this concurrent activity count as two separate eruptions (Kubanek et al., 2015)? It is difficult to answer these questions because of uncertainty in the timing of eruptive events, especially for those events observed solely in the geologic record, and it is therefore important to define eruptive events based on mapped relationships for long-term hazard assessments.

Exploring the sensitivity of the hazard assessment to the definition of an event requires additional assumptions about how to model event activity. Lava flow simulation depends on the volume and thickness of lava flows, variables we obtain from published data, as a proxy for eruption magnitude (Connor et al., 2012). Effusion rates are often highest during the initial phases of an eruption, which results in the maximum distal flow extent being reached early on (Bonny and Wright, 2017). Subsequent activity, marked by a lower effusion rate, is typically characterized by small length to width ratios (Kilburn and Lopes, 1988). Eruptions on the ESRP have likely had high initial effusion rates during the first phases of activity and then continued to build compound flow fields as effusion rate drops (Kuntz et al., 1992). We therefore simulate the initial phase of event eruptions using the same flow parameters as eruptions from single vents and assume effusion from a single event point, rather than distributing lava from a random number of vents and building compound flow fields. This is because effusion rate and volume, followed by eruption duration, are the main controls on flow length (Rowland et al., 2005).


Event Modeling

Our method for grouping vents into events employs similar clustering techniques to Runge et al. (2014), but departs from their use of expert elicitation in favor of spatio-temporal relationships identified from geologic mapping. We define a vent as a localized source of effusive activity. By contrast, we define an eruptive event as a vent or group of closely spaced vents erupted over a relatively short time interval. An event represents the complete record of activity related to ascent of a magma body or the emplacement of a series of dikes (e.g., Hell’s Half Acre, ESRP). Additionally, an event may also represent several subsequent eruptions from the same cone or fissure (e.g., COM). We define an event’s location as the mean of the coordinates of its near-neighbor vents (Fig. 1). While COM is compositionally more evolved than the majority of ESRP lavas, it is also spatially isolated from contemporaneously erupted non-COM sources, so events are not defined by geochemical variation. Uncertainty about the number of independent events in this study arises primarily because one-third of identified vents have no radiometric age determinations and their stratigraphic relationship to dated units is ambiguous.

A total of 506 surface vents have been mapped on the ESRP; 355 have an assigned age through 14C, K-Ar, Ar-Ar, or paleomagnetic dating methods, and 151 are undated (Kuntz et al., 1994, 2007; Anderson and Liszewski, 1997). We identify groups of vents from the 355 dated vents that may have formed as part of one larger event based on their temporal proximity to one another using a nearest neighbor clustering algorithm. We define these temporally congruous clusters in a way that captures natural breaks in the cumulative distribution of vent ages, which is controlled by the rate of volcanism and the resolution of methods used to date activity (Fig. DR1 in the GSA Data Repository1). The result of this grouping is a set of 52 clusters whose constituent vents were formed in close temporal proximity to one another (<1500 yr).

An elliptical template is positioned at the center of each temporal cluster to further identify relationships based on spatial proximity (Fig. 2). The ellipse is 20 km × 10 km with the long axis striking 330°. The dimensions and orientation were selected based on mapped ESRP dikes, non-eruptive fissures, and tension cracks (Kuntz et al., 1992, 1994, 2007), which reflect the plane normal to the regional least principal stress direction. We note this governs the emplacement geometries of propagating dikes, and not the overall spatial distribution trend of volcanic centers on the ESRP. If any vent within the cluster resides outside of the template, the cluster is broken into sub-clusters and templates are fit to the centers of these new clusters. The process repeats for all clusters until each vent resides within a template. The center of each of these templates is reported as the coordinate of each eruptive event. Vents without ages were also organized into events, independent of the dated events, using this template. The 355 vents with determined ages were grouped into 159 events, while the 151 undated vents were grouped into 97 events (Fig. 2) (Table DR1).

Spatial Density Estimation

Vent and event distributions are used to forecast the locations of potential volcanic eruptions on the ESRP. Burial of eruptive centers by lavas and sediment occurs non-uniformly on the ESRP, which biases spatial distribution models. This is particularly pronounced in basins due to a combination of non-uniform subsidence across the ESRP and burial that obscures the local eruptive history at an accelerated rate (Wetmore et al., 2009) (Fig. 1). We therefore include 32 buried eruptive centers, identified by Anderson and Liszewski (1997) and Wetmore et al., (2009), for both vent and event spatial density calculations to correct for some of this bias.

The spatial density of eruptive centers, the conditional probability of where a new vent will form, given that one forms somewhere on the ESRP, is estimated using a statistical model called nonparametric kernel density estimation (Connor and Connor, 2009; Bebbington and Cronin, 2011). We use a best-fit bivariate Gaussian kernel function with a directional smoothing bandwidth. The size, shape, and orientation of the kernel is determined by the locations of eruptive centers on the ESRP and not the regional alignment of dikes (Wetmore et al., 2009). The best-fit kernels for both the vents and events are elongate to the northeast, parallel to the overall trend of the ESRP (Fig. 3). Results show that areas of highest vent/event density correlate with the thickest total basalt distribution, suggesting that our modeled data effectively approximate spatial variations in long-term magma generation (Shervais et al., 2012).

Lava Flow Simulation

MOLASSES (MOdular LAva Simulation Software for Earth Science), a lava flow simulator modified from the LavaPL algorithm of Connor et al. (2012), distributes lava between cells based on rules that govern flow behavior (Kubanek et al., 2015). MOLASSES has been successfully benchmarked (Dietterich et al., 2017), performs well at recreating flow geometries similar to those found on the ESRP (Fig. DR2), and is sensitive to the geometries of lava flows, their thickness, area, and the underlying topography, rather than to the mechanics of lava flow emplacement. MOLASSES is useful for simulating the eventual footprint of a lava flow, but not its emplacement rate. Different types of lava flows result in different geometries and no single simulator is yet fully capable of modeling the complexities of all lava flow morphologies. We concentrate on the area inundated, recognizing these model limitations.

Inputs for MOLASSES include pulse volume, flow thickness, erupted volume, eruption location, and a digital elevation model (DEM) of the region. No spatial trend exists in the distribution of lava volumes or thicknesses across the ESRP, thus input parameters were stochastically sampled from probability distributions (Table DR2). Monte Carlo simulations onto a 90 m Shuttle Radar Topography Mission DEM (reset to the original topography for each flow) generated a range of conditional probabilities of site inundation of INL. The vent and event spatial density maps, along with the source locations that produced inundating flows, were used to identify areas of greatest hazard.

Of the 10,000 vent simulations, 3114 breached the INL border and 2024 initiated within its boundaries (Fig. 3). Additionally, 10,000 flows were simulated for event eruptions; 3209 of these flows partially inundated INL, with 2339 events initiating within its boundaries. Eruptive centers >30 km from the boundary of INL did not produce an inundating flow for either set of simulations. The probability of partial inundation of INL is ∼30%, given a future eruption. The conditional probability of lava inundation of INL, given an eruption in the region, is not particularly sensitive to event definition.

Recurrence Rate of Volcanism

The probability of lava flow inundation is made unconditional by accounting for the rate of volcanic activity. The recurrence interval between eruptions contributes to uncertainty in inundation probability calculations because it relies upon eruption catalog completeness and the accuracy of dating techniques. Several questions arise when selecting the appropriate data set for calculating the recurrence interval: does the eruption rate on the ESRP change with time? Is bias in eruption rate introduced through surface mapping and sampling due to burial of older lava flows by younger eruptions? Is additional bias introduced due to uncertainty in radiometric age determinations? Because the likely answer to these questions is ‘yes’, we must consider a variety of approaches in addressing how bias is introduced in the calculation of a recurrence interval.

The interval between eruptions is 2400 yr for mapped vents and is modeled at 4700 yr between events (Table DR3). An examination of the cumulative mapped vent count versus age suggests that recurrence rate was relatively constant from 500 ka through the beginning of the Holocene, after which COM initiated (Fig. DR1). It is likely that the estimated recurrence interval for activity older than 500 ka is biased due to burial by younger flows and sediments (Wetmore et al., 2009). We therefore take into account only the activity from 500 ka through the present for consideration in calculating the recurrence interval of volcanism on the ESRP (Fig. DR1). Using intervals of 1740 yr between eruptions for vents and 3800 yr for events, the annual probability of partial lava flow inundation of INL varies from 8.4 × 10−5 to 1.8 × 10−4. The annual probability of initiation of an eruption within the INL varies from 6.2 × 10−5 to 1.2 × 10−4 (Table 1).


Eruptive centers on the ESRP are generally well exposed, but variable subsidence and sedimentation obscure the most recent volcanism in some areas. Spatial density estimations that include buried eruptive centers in depocenters aids in removing some of the bias introduced by these processes (George et al., 2015). Like other volcanic fields, the geometries of eruptive centers on the ESRP vary from single vents to fissures and shields, creating a challenge for transforming map data of volcanic vents into a record of volcanic events. Yet, this transformation is essential for using geologic map data with simulation of lava flows and recurrence intervals of volcanism.

Our method provides a robust approach for addressing these issues, which are widespread in volcanic hazard assessments. The mapped geologic data are used to generate input parameters for the lava simulator, which is used to model the expected footprint of lava flows, rather than the dynamics of their emplacement. This output is used to determine the probabilities of lava flow inundation through Monte Carlo simulation. Results are coupled with recurrence intervals to calculate the annual unconditional probabilities of lava flow inundation and vent/event formation within the volcanic field.

At INL, relatively high conditional probability arises due to the position of the site in an area of spatially dense volcanic activity and its location in a topographic low, which tends to focus lava flows from vents outside the INL boundaries onto the property. We estimate this conditional probability to be ∼30% for the entire site, which exceeds International Atomic Energy Agency guidelines for nuclear facilities (IAEA, 2016). Volcanic risk to individual facilities could be estimated by using higher-resolution DEMs and site-specific engineering data.


Funding provided by National Science Foundation grant ACI 1339768 SI2-SSI: Collaborative Research: Building Sustainable Tools and Collaboration for Volcanic and Related Hazards. Publication funding generously provided by the USF libraries and School of Geosciences. We thank Mark Quigley (editor) and three anonymous reviewers, whose comments improved this paper.

1GSA Data Repository item 2018326, Table DR1 (lava flow model variables), Table DR2 (ESRP event clusters), Table DR3 (ESRP recurrence intervals), Figure DR1 (plot of cumulative vent distribution over time), and Figure DR2 (comparison of mapped and simulated lava flows), is available online at http://www.geosociety.org/datarepository/2018/ or on request from editing@geosociety.org.
Gold Open Access: This paper is published under the terms of the CC-BY license.