- Methodology
- Open Access
Probabilistic approach to modeling lava flow inundation: a lava flow hazard assessment for a nuclear facility in Armenia
- Laura J Connor^{1}Email author,
- Charles B Connor^{1},
- Khachatur Meliksetian^{2} and
- Ivan Savov^{3}
https://doi.org/10.1186/2191-5040-1-3
© Connor et al; licensee Springer. 2012
- Received: 20 June 2011
- Accepted: 25 January 2012
- Published: 25 January 2012
Abstract
Probabilistic modeling of lava flow hazard is a two-stage process. The first step is an estimation of the possible locations of future eruptive vents followed by an estimation of probable areas of inundation by lava flows issuing from these vents. We present a methodology using this two-stage approach to estimate lava flow hazard at a nuclear power plant site near Aragats, a Quaternary volcano in Armenia.
Keywords
- lava flow simulation
- modeling code
- probabilistic hazard assessment
- spatial density
- Monte Carlo method
- Armenia
Background
Volcanic hazard assessments are often conducted for specific sites, such as nuclear facilities, dams, ports and similar critical facilities that must be located in areas of very low geologic risk (Volentik et al 2009; >Connor et al 2009). These hazard assessments consider the hazard and risk posed by specific volcanic phenomena, such as lava flows, tephra fallout, or pyroclastic density currents (IAEA 2011; Hill et al 2009). Although site hazards could be considered in terms of the cumulative effects of these various volcanic phenomena, a better approach is to assess the hazard and risk of each phenomenon separately, as they have varying characteristics and impacts. Here, we develop a methodology for site-specific hazard assessment for lava flows. Lava flows are considered to be beyond the design basis of nuclear facilities, meaning that the potential for the occurrence of lava flows above some level of acceptable likelihood would exclude the site from development of nuclear facilities because safe control or shutdown of the facility under circumstances of lava flow inundation cannot be assured (IAEA 2011).
Spatial density estimation
Site-specific lava flow hazard assessments require that the hazard of lava inundation be estimated long before lava begins to erupt from any specific vent. In many eruptions, lavas erupt from newly formed vents, hence, the potential spatial distribution of new vents must be estimated as part of the analysis. This is particularly important because the topography around volcanoes is often complex and characterized by steep slopes. Small variations in vent location may cause lava to flow in a completely different direction down the flanks of the volcano. Thus, probabilistic models of lava flow inundation are quite sensitive to models of vent location. Furthermore, many volcanic systems are distributed. Examples include monogenetic volcanic fields (e.g. the Michoacán-Guanajuato volcanic field, Mexico), distributed composite volcanoes which lack a central crater (e.g. Kirishima volcano, Japan), and volcanoes with significant flank activity (e.g. Mt. Etna, Italy). Spatial density estimates are also needed to forecast potential vent locations within such distributed volcanic systems (Cappello et al 2011).
In addition, loci of activity may wax and wane with time, such that past vent patterns may not accurately forecast future vent locations (Condit and Connor 1996). Thus, it is important to determine if temporal patterns are present in the distribution of past events, so that an appropriate time interval can be selected for the analysis (i.e., use only those vents that represent likely future patterns of activity, not older vents that may represent past patterns).
Kernel density estimation is a non-parametric method for estimating the spatial density of future volcanic events based on the the locations of past volcanic events (Connor and Connor 2009; Kiyosugi et al 2010; Bebbington and Cronin 2010). Two important parts of the spatial density estimate are the kernel function and its bandwidth, or smoothing parameter. The kernel function is a probability density function that defines the probability of future vent formation at locations within a region of interest. The kernel function can be any positive function that integrates to one. Spatial density estimates using kernel functions are explicitly data driven. A basic advantage of this approach is that the spatial density estimate will be consistent with known data, that is, the spatial distribution of past volcanic events. A potential disadvantage of these kernel functions is that they are not inherently sensitive to geologic boundaries. If a geologic boundary is known it is possible to modify the density estimate with data derived from field observations and mapping. Connor et al (2000) and Martin et al (2004) discuss various methods of weighting density estimates in light of geological or geophysical information, in a manner similar to Ward (1994). A difficulty with such weighting is the subjectivity involved in recasting geologic observations as density functions.
Equation 1 is a simplification of this more general case, whereby the amount of smoothing by the bandwidth, h, varies consistently in both the N-S and E-W directions. The bandwidth, H, on the other hand, is a 2 × 2 element matrix that specifies two distinct smoothing patterns, one in a N-S trending direction and another in an E-W trending direction. This bandwidth matrix is both positive and definite, important because the matrix must have a square root. |H| is the determinant of this matrix and H^{-1/2} is the inverse of its square root. x is a 1 × 2 distance matrix (i.e. the x-distance and y-distance from s to an event), b is the cross product of x and H^{-1/2}, and b^{ T } is its transform. The resulting spatial density at each point location, s, is usually distributed on a grid that is large enough to cover the entire region of interest. Bandwidth selection is a key feature of kernel density estimation (Stock and Smith 2002; Connor et al 2000; Molina et al 2001; Abrahamson 2006; Jaquet et al 2008; Connor and Connor 2009), and is particularly relevant to lava flow hazard studies. Bandwidths that are narrow focus density near the locations of past events. Conversely, a large bandwidth may over-smooth the density estimate, resulting in unreasonably low density estimates near clusters of past events, and overestimate density far from past events. This dependence on bandwidth can create ambiguity in the interpretation of spatial density if bandwidths are arbitrarily selected. A further difficulty with elliptical kernels is that all elements of the bandwidth matrix must be estimated, that is the magnitude and direction of smoothing in two directions. Several methods have been developed for estimating an optimal bandwidth matrix based on the locations of the event data (Wand and Jones 1995), and have been summarized by Duong (2007). Here we utilize a modified asymptotic mean integrated squared error (AMISE) method, developed by Duong and Hazelton (2003), called the SAMSE pilot bandwidth selector, to optimally estimate the smoothing bandwidth for our Gaussian kernel function. These bandwidth estimators are found in the freely available R Statistical Package (Hornik 2009; Duong 2007). Bivariate bandwidth selectors like the SAMSE method are extremely useful because, although they are mathematically complex, they find optimal bandwidths using the actual data locations, removing subjectivity from the process. The bandwidth selectors used in this hazard assessment provide global estimates of density, in the sense that one bandwidth or bandwidth matrix is used to describe variation across the entire region.
Given that spatial density estimates are based on the distribution of past volcanic events, existing volcanic vents within a region and time period of interest first need to be identified and located. This compilation is then used as the basis for estimating the probability of the opening of new vents within a region. Our lava flow hazard assessment method is concerned with the likelihood of the opening of new vents that erupt lava flows. Such vents may form when magma first reaches the surface, forming a new volcano, or may form during an extended episode of activity, whereby multiple vents may form while an eruptive episode continues over some period of time, generally months to years (Luhr and Simkin 1993), and the locus of activity shifts as new dikes are injected into the shallowest part of the crust. Therefore, for the purposes of this study, an event is defined as the opening of a new vent at a new location during a new episode of volcanic activity. Multiple vents formed during a single episode of volcanism are not simulated.
Numerical Simulation of Lava flows
On land, a lava flow is a dynamic outpouring of molten rock that occurs during an effusive volcanic eruption when hot, volatile-poor, relatively degassed magma reaches the surface (Kilburn and Luongo 1993). These lava flows are massive volcanic phenomena that inundate areas at high temperature (> 800°C), destroying structures, even whole towns, by entombing them within meters of rock. The highly destructive nature of lava flows demands particular attention when critical facilities are located within their potential reach.
The area inundated by lava flows depends on the eruption rate, the total volume erupted, magma rheological properties, which in turn are a function of composition and temperature, and the slope of the final topographic surface (Dragoni and Tallarico 1994; Griffiths 2000; Costa and Macedonio 2005). Previous studies have modeled the physics of lava flows using the Navier-Stokes equations and simplified equations of state (Dragoni 1989; Del Negro et al 2005; Miyamoto and Sasaki 1997). Other studies have concentrated on characterizing the geometry of lava flows, and studying their development during effusive volcanic eruptions (Walker 1973; Kilburn and Lopes 1988; Stasiuk and Jaupart 1997; Harris and Rowland 2009). These morphological studies are mirrored by models that concentrate on the areal extent of lava flows, rather than their flow dynamics. These models generally abstract the highly complex rheological properties of lava flows using geometric terms and/or simplified cooling models (Barca et al 1994; Wadge et al 1994; Harris and Rowland 2001; Rowland et al 2005).
A new lava flow simulation code, written in PERL, was created to assess the potential for site inundation by lava flows, similar, in principle, to areal-extent models. This lava flow simulation tool is used to assess the probability of site inundation rather than attempting to model the complex real-time physical properties of lava flows. Since the primary physical information available for lava flows is their thickness, area, length and volume, this model is guided by these measurable parameters and not directly concerned with lava flow rates, their fluid-dynamic properties, or their chemical makeup and composition. The purpose of the model is to determine the conditional probability that flow inundation of a site will occur, given an effusive eruption at a particular location estimated using the spatial density model discussed previously.
A total volume of lava to be erupted is set at the start of each model run. The model assumes that each cell inundated by lava retains or accumulates a residual amount of lava. The residual must be retained in a cell before that cell will pass any lava to adjacent cells. This residual corresponds to the modal thickness of the lava flow. Lava may accumulate in any cell to amounts greater than this residual value if the topography allows pooling of lava. As flow thickness varies between lava flows, the residual value chosen for the flow model also varies from simulation to simulation. Here, our term residual corresponds to the term adherence, used in codes developed by Wadge et al (1994) and Barca et al (1994). In our case, residual lava does not depend on temperature or underlying topography, but rather, is used to maintain a modal lava flow thickness. Lava flow thicknesses, measured within the site area, are fit to a statistical distribution which is sampled stochastically in order to choose a residual (i.e. modal thickness) value for each realization. Lava flow simulation requires a digital elevation model (DEM) of the region of interest. One source of topographic DEM data is the Shuttle RADAR Topography Mission (SRTM) database. The 90-meter grid spacing of SRTM data limits the resolution of the lava flow. Topographic details smaller than 90 m can influence flow path, but these cannot be accounted for using a 90-m DEM. A more detailed DEM could provide enhanced flow detail, but a decrease in DEM grid spacing increases the total number of grid cells, thus increasing computation time as the flow has to pass through an increasing number of grid cells. A balance needs to be maintained between capturing important flow detail over the topography and limiting the overall time required to calculate the full extent of the flow. Critical considerations for grid spacing are the topography of the site area and the volumes and flow rates of local lava flows. Lava flows erupted at high rate or high viscosity would quickly overwhelm surrounding topography, so in these cases a coarse 90-m DEM may be sufficient for flow modeling. For low flow rates or low viscosities, lava flows would meander around smaller topographic features which would be unresolved in a coarse 90-m DEM. Therefore, in these cases a higher resolution DEM would be necessary to achieve credible model results. In our study, a 90-m DEM was considered adequate due to the unavailability of information regarding lava flow rates in the area and assumed higher flow rates based on flow geometries measured in the field. Also, the boundaries of the plateau on which the ANPP site is located was determined to be adequately resolved by a 90-m DEM.
Lava flow hazard at the Armenian nuclear power plant site
Lava flows are a common feature of the Armenian landscape. Some mapped flows are highlighted in Figure 2. A group of 18 volcanic centers comprise an area known as the Shamiram Plateau (this area is located within the red box in Figure 1). The Armenian nuclear power plant (ANPP) site lies within this comparatively dense volcanic cluster at the southern margin of the Shamiram Plateau. Our lava flow hazard assessment is designed to assess the conditional probability that lava flows reach the boundary of the site area, given an effusive eruption on the Shamiram Plateau. In addition, large-volume lava flows are found on the flanks of Aragats volcano, a 70-km-diameter basalt-trachyandesite to trachydacite volcano located immediately north of the Shamiram Plateau.
The mapped lava flows on the Shamiram Plateau can be divided into two age groups, pre-ignimbrite lava flows that range in age from approximately 0.91-1.1 Ma, and post-ignimbrite lava flows that cover the ignimbrites of Aragats volcano. The youngest features of Aragats Volcano are large volume lava flows from two cinder cones, Tirinkatar (0.45 Ma) and Ashtarak (0.53 Ma). All of these age determinations are based on K-Ar dating by Chernyshev et al (2002). The youngest small-volume lava flows of the Shamiram Plateau are the Dashtakar group of cinder cones, based on borehole evidence indicating that the Dashtakar flows overlay one of these ignimbrites of Aragats.
Lava flows of the Shamiram Plateau are typical of monogenetic fields, being of comparatively low volume, generally < 0.03 km^{3}, and short total length, generally < 5 km. Based on logging data from four boreholes and including the entire area of the Shamiram Plateau and estimated thickness of the lava pile, the total volume of lava flows making up the plateau is ~11- 24 km^{3}. Given these values, hundreds of individual lava flows comprise the entire plateau. Thus, there is a possibility that lava flows will inundate the site in the future, associated with the eruption of monogenetic volcanoes on the Shamiram Plateau, should such eruptions occur.
Size estimates of lava flows
Volcano (source) | Area (km^{2}) | Thickness (m) | Volume (km^{3}) | Length (km) | Composition |
---|---|---|---|---|---|
Arich | 16.3 | 8 | 0.130 | 9.48 | TB^{1}, BTA^{1} |
Atomakhumb | 3.9 | 6 | 0.023 | 3.43 | BA^{1}, BTA |
Barcradir(Bartsradir) | 32.9 | 9 | 0.296 | 12.10 | TB, BTA |
Bazmaberd | 13.1 | 14 | 0.184 | 6.34 | BA, BTA |
Blrashark | 1.6 | 6 | 0.010 | 2.49 | TA^{1}, TD^{1} |
Blrashark | 2.5 | 7 | 0.018 | 3.13 | TA, TD |
Bolorsar | 2.2 | 6 | 0.013 | 2.72 | BTA, TA |
Dashtakar | 2.1 | 10 | 0.021 | 4.44 | BA, BTA |
Dashtakar | 1.6 | 6 | 0.009 | 3.66 | BA, BTA |
Karmratar | 0.7 | 4 | 0.003 | 3.61 | TA |
Mets Mantash | 8.9 | 9 | 0.080 | 8.47 | TB, BTA |
Shamiram | 1.0 | 4 | 0.004 | 1.41 | TA |
Siserasar | 0.8 | 11 | 0.009 | 1.72 | TA |
Tirinkatar-2 | 13.3 | 4 | 0.053 | 6.54 | BTA, BA |
Topqar(Topkar) | 2.9 | 9 | 0.026 | 3.07 | BTA, TA |
Ashtarak | 84 | 6 | 0.50 | 26.50 | BA, BTA |
Irind | 66 | 55 | 3.65 | 20.53 | Dacite |
Paros | 109 | 8 | 0.87 | 33.36 | TB, BTA |
Tirinkatar-1 | 75 | 7 | 0.53 | 26.36 | BTA, BA |
Pokr Bogutlu | 165 | 110 | 18.18 | 27.92 | TD |
(Cakhkasar) |
Volume-limited flows occur when small batches of magma reach the surface and erupt for a brief period of time, forming lava flows associated with individual monogenetic centers. These eruptions often occur in pulses and erupting vents may migrate a short distance, generally < 1 km, during the eruption. Each pulse of activity in the formation of the monogenetic center may produce a new individual lava flow, hence, constructing a flow field over time. The longest lava flows in these fields are generally those associated with the early stages of the eruption, when eruption rates are greatest (Kilburn and Lopes, 1988). Within the Shamiram Plateau area, individual monogenetic centers have one (e.g. Shamiram volcano) to many (e.g. Blrashark volcano) individual lava flows.
Longer lava flows are also found on Aragats volcano, especially higher on its flanks (Table 1). These summit lavas comprise a thick sequence of trachyandesites and trachydacites having a total volume > 500 km^{3}. The most recent lava flows from the flanks of Aragats include Tirinkatar, which is separated into two individual trachybasalt flows Tirinkatar-1 and Tirinkatar-2, and the Ashtarak lava flow. Tirinkatar-1 and Ashtarak each have volumes ~0.5 km^{3}. The largest volume flank lava flows are part of the trachydacitic Cakhkasar lava flow of Pokr Bogutlu volcano, with a total volume ~18 km^{3}, on the same order as the largest historical eruptions of lava flows worldwide (Thordarson and Self 1993). These larger volume lava flows are effusion rate-limited, since the length of the lava flow is controlled by the effusion rate at the vent. The lengths of the Ashtarak and Tirinkatar-1 lava flows exceed 20 km. Based on comparison with observed historical eruptions, their effusion rates were likely on the order of 100 m^{3} s^{-1} (Walker, 1973; Malin 1980; Kilburn and Lopes, 1988; Harris and Rowland, 2009). Thus, while volume-limited flows erupt on the Shamiram Plateau in the immediate vicinity of the site, effusion rate-limited flows erupt at higher elevations on the flanks of Aragats volcano. While it is conceivable that these larger volume flows may reach the site because of their great potential length, this event is less likely because their occurrence is so infrequent. Another deterrent is the fact that the Shamiram plateau acts as a topographic barrier to these longer, larger flows reaching the ANPP site. Each class of lava flows, smaller volume-limited flows and larger effusion rate-limited flows, is considered separately when assessing lava flow hazard at the ANPP site.
Results and Discussion
Using spatial density estimation
Locating the source region of erupting lava is critical in determining the area inundated by a lava flow. Probable source regions are estimated using a spatial density model, which in turn depends on a geological map identifying the locations of past eruptive vents. In this context, volcanic vents are defined as the approximate locations where magma has or may have reached the surface and erupted in the past. A primary difficulty in using a data set of the distribution of volcanic vents is determination of independence of events. In statistical parlance, independent events are drawn from the same statistical distribution, but the occurrence of one event does not influence the probability of occurrence of another event. We are interested in constructing a spatial density model only using independent events. Unfortunately, it is difficult to determine from mapping and stratigraphic analysis if vents formed during the same eruptive episode or occurred as independent events during different volcanic eruptions. Some of these are easily recognized (e.g. boccas that are located adjacent to scoria cones). In other cases, it is uncertain if individual volcanoes should be considered to be independent events, or were in reality part of the same event. Because of this uncertainty, alternative data sets are useful when estimating the spatial density. Here, we use one data set to maximize the potential number of volcanic events: all mapped vents are included in the data set as independent events. An alternative data set could consider volcanic events to be comprised of groups of volcanic vents that are closely spaced and not easily distinguished stratigraphically.
In order to apply the spatial density estimate, it is assumed that 18 mapped volcanic centers represent the potential distribution of future volcanic vents on the Shamiram Plateau. Some older vents are no doubt buried by subsequent volcanic activity. It is also possible that older vents are buried in sediment of the Yerevan basin, south of the ANPP site.
Volcanic vents mapped on the Shamiram Plateau
Easting | Northing |
---|---|
425507 | 4449732 |
425649 | 4449144 |
425992 | 4449400 |
425053 | 4449362 |
428682 | 4452894 |
429363 | 4452946 |
429504 | 4452711 |
429931 | 4452251 |
427322 | 4449676 |
427383 | 4449840 |
427835 | 4450008 |
428332 | 4444255 |
427386 | 4454344 |
427538 | 4453062 |
430618 | 4442102 |
427623 | 4452343 |
426857 | 4451520 |
425285 | 4454652 |
A grid-based flow regime
The SRTM database from CGIAR-CSI (the Consultative Group on International Agricultural Research-Consortium for Spatial Information) is used as a model of topographic variation on the Shamiram Plateau and adjacent areas. This consortium (Jarvis et al, 2008) has improved the quality of SRTM digital topographic data by further processing version 2 (released by NASA in 2005) using hole-filling algorithms and auxiliary DEMs to fill voids and provide continuous topographical surfaces. For the lava flow simulation, these data are converted to a UTM Zone 38 N projection, using the USGS program, PROJ4, and re-sampled at a 100 × 100 m grid spacing, using the mapping program GMT. In the model, lava is distributed from one 100 m^{2} grid cell to its adjacent grid cells.
The region that was chosen for the lava flow model is identified in Figure 1 (red-dashed box). Within this area a new vent location is randomly selected based on a spatial density model of 18 events clustered within and around the Shamiram Plateau (Figure 4). The model simulates a flow of lava from this new vent location onto the surrounding topography. The total volume of lava to be erupted is specified at the onset of a model run. Lava is added incrementally to the DEM surface at the vent location until the total specified lava flow volume is reached. At each iteration, 10^{5} m^{3} is added to the grid cell at the location of the vent (source) and is distributed over adjoining grid cells. Given that a grid cell is 100 m^{2}, this corresponds to adding a total depth of 10 m to the vent cell at each iteration.
The lava flow simulation is not intended to mimic the fluid-dynamics of lava flows, so these iterations are only loosely associated with time steps. For example, volume-limited lava flows of the Shamiram Plateau are generally < 5 km in length, with volumes on the order of 0.3 - 2.3 × 10^{-2} km^{3}. These volumes and lengths agree well with lavas from compilations by Malin (1980) and Pinkerton and Wilson (1994). For such lava flows, effusion rates of 10 - 100 m^{3} s^{-1} are expected (Harris and Rowland, 2009). Using these empirical relations, an iteration adding a volume of 10^{5} m^{3} of lava corresponds to an elapsed time of 10^{3}- 10^{4}s. Lava is distributed to adjacent cells only at each iteration, so this effusion rate corresponds to flow-front velocity on the order of 0.01 - 0.1 ms ^{ - }^{1}, in reasonable agreement with observations of volume-limited flow-front velocities.
Parameter estimation for Monte Carlo simulation
Lava flow simulation input parameters
Parameter | Range | Notes |
---|---|---|
ANPP site boundary | Boundaries used in analysis | |
East (km) | 428.2 | |
West (km) | 426.0 | |
North (km) | 4449.0 | |
South (km) | 4447.0 | |
Lava thickness (m) | 4-15 | Truncated normal distribution; |
Mean = 7.0 m | ||
Standard Dev. = 3.0 m | ||
Lava flow volume (m^{3}) | 10^{6}-10^{9} | Truncated normal distribution; |
(log)Mean = 7.2 | ||
(log)Standard Dev. = 0.5 | ||
Iteration volume | 10^{5} | Lava volume added at source |
vent in each iteration | ||
Number of simulations | 10 000 |
Configuration file for lava flow simulation of vents on the Shamiram Plateau
Parameter = Value | Explanation |
---|---|
Inputs | |
DEM_SOUTH = 4440 | N, S, E, W |
DEM_NORTH = 4470 | boundaries |
DEM_EAST = 440 | of the DEM |
DEM_WEST = 410 | |
DEM_SPACING = 0.1 | DEM grid spacing (km) |
DEM_FILE = file | (ASCII format) rows of elevation values (masl) |
RESIDUAL_AV = 8.0 | Lava thickness (m): Average |
RESIDUAL_SD2 = 1.0 | Standard Deviation |
(higher value=higher lava viscosity) | |
ERUPTED_LAVA = 1e5 | Volume of lava distributed |
per iteration or pulse (m^{3}) | |
TOTAL_LAVA_AV = 1e7 | Lava volume (m^{3}): Average |
TOTAL_LAVA_SD2 = 0.5 | Standard Deviation |
FLOWS = 1 | Number of lava flows to simulate per run |
RUNS = 10 000 | Number of lava flow runs (for statistical analysis) |
AOI_WEST = 426.0 | Area of interest |
AOI_EAST = 428.2 | |
AOI_SOUTH = 4447.8 | |
AOI_NORTH = 4449.0 | |
SPATIAL_DENSITY_FILE = file | X Y Z format, grid of spatial density |
values for the potential of volcanism | |
SPATIAL_DENSITY_SPACING=.1 | spacing of spatial density grid (km) |
In reality, more than one lava flow may erupt during the course of formation and development of a single monogenetic volcano. However, the first lava flow to form during this eruption will tend to have the longest length and greatest potential to inundate the ANPP site. Experiments were conducted to simulate the formation of multiple (up to 10) lava flows from a single vent, or group of closely spaced vents. It was determined that the later lava flows tend to broaden the flow field, but not lengthen it. This result is in agreement with observations of lava flow field development on Mt. Etna (Kilburn and Lopes, 1988). For the ANPP site, the conditional probability of site inundation was sensitive to lava flow length, but insensitive to broadening of the lava flow field. Therefore, only one lava flow was simulated per eruptive vent. Nevertheless, for some sites the potential for broadening the area of inundation by successive flows may be an important factor.
Simulation results
A total of 10 000 simulations were executed in order to estimate the probability of lava flow inundation resulting from the formation of new monogenetic vents on the Shamiram Plateau. Out of 10 000 events, 2485 of the simulated flows crossed the perimeter of the site, or 24.9% percent of the total number of simulations.
27 Mapped vents on the flanks of Aragats Volcano
Easting | Northing |
---|---|
430920 | 4485826 |
422295 | 4488512 |
414366 | 4498480 |
439898 | 4478024 |
440441 | 4476970 |
425896 | 4491003 |
421407 | 4471589 |
418534 | 4469462 |
408119 | 4495051 |
408990 | 4481638 |
414068 | 4471495 |
427253 | 4483296 |
424558 | 4482259 |
423136 | 4480327 |
411159 | 4469329 |
423682 | 4494414 |
405800 | 4477396 |
406683 | 4476948 |
418530 | 4494870 |
424111 | 4495248 |
408363 | 4492635 |
415964 | 4497175 |
422344 | 4491454 |
428042 | 4474090 |
428225 | 4474806 |
424775 | 4492714 |
399806 | 4491891 |
Configuration file for simulation of lava flows from flank vents
Parameter = Value | Explanation |
---|---|
Inputs | |
DEM_SOUTH = 4441 | N, S, E, W |
DEM_NORTH = 4482 | boundaries |
DEM_WEST = 408 | of the DEM |
DEM_EAST = 448 | |
DEM_SPACING = 0.1 | DEM grid spacing (km) |
DEM_FILE = file | (ASCII format) rows of elevation values (masl) |
MIN_RESIDUAL = 1 | Map to observed flow thicknesses (m): |
MAX_RESIDUAL = 4 | 1 = 6, 2 = 7, 3 = 8, 4 = 9, |
(lower value = lower lava viscosity) | |
ERUPTED_LAVA = 1e6 | Volume of lava distributed |
per iteration or pulse (m^{3}) | |
MIN_TOTAL_LAVA = 1 | Map to observed flow volumes (km^{3}): |
MAX_TOTAL_LAVA = 3 | 1 = 5 × 10^{8} 2 = 5.3 × 10^{8}, 3 = 8.7 × 10^{8} |
FLOWS = 1 | Number of lava flows to simulate per run |
RUNS = 1000 | Number of lava flow runs |
AOI_WEST = 426.0 | Area of interest |
AOI_EAST = 428.2 | |
AOI_SOUTH = 4447.8 | |
AOI_NORTH = 4449.0 | |
SPATIAL_DENSITY_FILE = file | X Y Z format, grid of spatial density |
values for the potential of volcanism | |
SPATIAL_DENSITY_SPACING=.1 | spacing of spatial density grid (km) |
Conclusions
We demonstrate a methodology for site-specific lava flow hazard assessment. This two-stage process uses a two-dimensional elliptical Gaussian kernel function to estimate spatial density. The SAMSE method, a modified asymptotic mean squared error approach, uses the distribution of known eruptive vents to optimally determine a smoothing bandwidth for the Gaussian kernel function. Potential vent locations (N = 10 000) are stochastically sampled from the resulting spatial density probability map. For each randomly sampled vent location, a lava flow inundation model is executed. Lava flow input parameters (volume and modal thickness) are determined from distributions fit to field observations of the low viscosity trachybasalt to trachydacite lava flows of the area. The areas and flow extents (a quantitative measure of lava flow length) of these simulated lava flows compare reasonably with those of mapped lava flows. This approach yields a conditional probability of lava flow inundation, given the opening of a new vent, and provides a map of vent locations leading to site inundation.
Lava flow hazards exist at the ANPP site because potential eruptions on the Shamiram Plateau may produce lava flows that inundate the site. This Monte Carlo analysis has shown that, given the number of relatively small-volume lava flows occurring on the Shamiram Plateau, approximately 25% of all eruptions, resulting from the formation of a new vent, might also produce lava flows that inundate the ANPP site. Although very long and voluminous lava flows occur in the Aragats volcanic system, this analysis demonstrates that these types of flows do not present a credible hazard for the site, as the topography of the Shamiram Plateau would divert such potential flows away from the site area.
An integrated hazard assessment also depends on the estimation of the recurrence rate of effusive volcanism. Assuming a recurrence rate of effusive eruptions on the Shamiram Plateau of 4.1 × 10^{-7} yr^{-1} and 3.5 × 10^{-6} yr^{-1}, based on currently available radiometric age determinations (Chernyshev et al, 2002), the annual probability of site inundation by renewed effusive volcanism on the Shamiram Plateau is approximately 1.0 × 10^{-7} to 8.8 × 10 ^{ - }^{7}.
Methods
Spatial density analysis
The lava flow hazard assessment begins with a spatial density analysis involving the locations of 18 volcanic events located on the Shamiram plateau. This analysis will help determine the most likely locations of future volcanic events which will then become the source locations for possible lava flows. These events are listed in Table 2. Using these 18 events an optimal bandwidth is determined using the SAMSE method in the 'ks' package within the statistical program, 'R'. The required 'R' commands are the following:
library (ks)
vents18 <- read.table ("events_zoom.wgs84.z38.utm")
show (vents18)
bw_samse_18vents <- Hpi(x = vents18, nstage = 2, pilot="samse", pre="sphere")
show (bw_samse_18vents)
where 'ks' is the name of the 'R' package needed to perform the analysis, vent18 is a local 'R' variable holding vent locations, vents _18 _wgs84.z38.utm is the input text file containing the vent locations (easting and northing separated by a space), bw _samse _18vents is a local 'R' variable holding the output from the 'Hpi' routine, the bandwidth matrix in meters:
[,1] [,2]
[1,] 844328.34 -13235.75
[2,] -13235.75 2113393.17
Spatial density analysis is accomplished using a PERL script (see Additional file 1). Parameters for the script are inserted directly at the top of the script as shown in the following code section:
#############################################################
# INPUT SECTION: These variables can be adjusted by the user
############################################################
## This is the complete set of events:
# events_zoom.wgs84.z38.utm: N = 18 <425053/430618> <4442102/4454652>
my $west = 420000;
my $east = 436000;
my $south = 4439000;
my $north = 4463000;
my $Grid_spacing = 100;
# The bandwidth matrix via SAMSE 2-stage pre-transformation 'sphering'
# units = square meters, for 18 events near ANPP
# [,1] [,2]
# [1,] 844328.34 -13235.75
# [2,] -13235.75 2113393.17
# units = square kilometers
my $H = pdl [
[.84432834, -.01323575],
[-.01323575, 2.11339317]
];
# The input file of event locations
my $in = "events_zoom.wgs84.z38.utm";
# The output file for the spatial intensity grid
my $out1 = "spatial_density_samse_events_zoom.wgs84.z38.utm.2";
where, $north, $south, $east, $west are the map boundaries in UTM meters, $Grid _spacing is the map grid spacing, units in km, $H is the kernel bandwidth, units converted to km^{2}, $in is the name of the input file of volcanic event locations (ASCII format: easting northing), and $out1 is the name of the output file of the spatial density grid (ASCII format: easting northing density).
$H is a matrix and its structure in the script is controlled by the PERL package 'pdl'. The 4 values for the matrix are derived from the output of the 'Hpi' routine (as noted above). To run the script from the command line type:
perl gausXY.pl
where 'gausXY.pl' is the name of the script. All parameters are inserted directly at the top of the script as indicated above.
A second PERL script drives the lava flow simulation (see Additional file 2 and Additional file 3). The inputs for this script are contained in a configuration file. To run the code from the command line type:
perl lava_flow.pl lavaflow.conf 0
where 'lava_flow.pl' is the name of the script, 'lavaflow.conf' is the name of the configuration file, and 0 is the starting run number. Each run of the script simulates one complete lava flow simulation. The total number of simulated lava flows is set in the configuration file. The configuration file parameters are listed in Tables 4 and 6.
The PERL lava flow simulation script produces 3 output files:
lava_flow_stats.# This file is a compilation of all simulated lava flows (where '#' refers to the initial run number). This text file contains 6 columns:
Easting (units = km)
Northing (units = km)
Hit (1 = hit; 0 = miss)
TL (units = cubic meters)
Residual (units = meters)
Total (units = cubic meters)
where Easting and Northing refer to the location of the erupting vent, Hit is either 1 or 0, where 1 means that the lava flow penetrated the area of interest (i.e. the boundary of the site) and 0 indicates that it did not, TL is the total volume of erupted lava, Residual refers to a flow's modal thickness, and Total is also the total volume erupted, but calculated in a different way. The total number of lava flow simulations are recorded.
flow.#.utm This file records the grid location and thickness of lava in each inundated cell (where '#' refers to an individual run number). This text file contains 3 columns: X Y thickness, where X Y refers to the inundated grid cell, and thickness refers to the thickness (m) of lava in that cell. This file is used to calculate the length and area of each simulated lava flow.
vents.utm This text file records the vent location of each lava flow simulation. The file contains two columns: Easting Northing.
Declarations
Acknowledgements
The authors gratefully acknowledge the logistical and technical support of Staff at the Institute of Geological Sciences of Armenian National Academy of Sciences. Discussions with Arkadi Karakhanian regarding Armenian geology and field mapping greatly enhanced the authors' overall understanding of the geological setting of Armenia. Reviews of early results of this study by Britt Hill, Willy Aspinall, and Antonio Godoy, all representing the International Atomic Energy Agency, led to improvements in the methods presented here. This research was partially supported by a grant from the US National Science Foundation (DRL 0940839). Reviews by Britt Hill and Antonio Costa improved the manuscript.
Authors’ Affiliations
References
- Abrahamson N: Seismic hazard assessment: Problems with current practice and future developments. First European Conference on Earthquake Engineering and Seismology, Geneva, Switzerland 2006.Google Scholar
- Barca D, Crisci GM, Gregorio SD, Nicoletta F: Cellular automata for simulating lava flows: A method and examples of the Etnean eruptions. Transport Theory and Statistical Physics 1994, 23: 195–232. 10.1080/00411459408203862View ArticleGoogle Scholar
- Bebbington MS, Cronin SJ: Spatio-temporal hazard estimation in the Auckland volcanic field, New Zealand, with a new event-order model. Bulletin of Volcanology 2010,73(1):55–72.View ArticleGoogle Scholar
- Cappello A, Vicari A, Del Negro C: Assessment and modeling of lava flow hazard on Mt. Etna volcano. Bollettino di Geofisica Teorica ed Applicata 2011,52(2):10–20.Google Scholar
- Chernyshev IV, Lebedev VA, Arakelyants MM, Jrbashyan R, Ghukasyan Y: Geochronology of the Aragats volcanic center, Armenia: Evidence from K-Ar dating. Doklady Earth Sciences 2002,384(4):393–398. (in Russian) (in Russian)Google Scholar
- Condit CD, Connor CB: Recurrence rate of basaltic volcanism in volcanic fields: An example from the Springerville Volcanic Field, AZ, USA. Geological Society of America Bulletin 1996, 108: 1225–1241. 10.1130/0016-7606(1996)108<1225:RROVIB>2.3.CO;2View ArticleGoogle Scholar
- Connor CB, Connor LJ: Estimating spatial density with kernel methods. In Volcanic and Tectonic Hazard Assessment for Nuclear Facilities. Edited by: Connor C, Chapman N, Connor L. Cambridge University Press; 2009:331–343.View ArticleGoogle Scholar
- Connor CB, Hill BE: Three nonhomogeneous Poisson models for the probability of basaltic volcanism: Application to the Yucca Mountain region. Journal of Geophysical Research 1995, 100: 12 107–10 125.Google Scholar
- Connor CB, Stamatakos JA, Ferrill DA, Hill BE, Ofoegbu GI, Conway FM, Sagar B, Trapp J: Geologic factors controlling patterns of small-volume basaltic volcanism: Application to a volcanic hazards assessment at Yucca Mountain, Nevada. Journal of Geophysical Research 2000,105(1):417–432. 10.1029/1999JB900353View ArticleGoogle Scholar
- Connor CB, Sparks RSJ, Díez M, Volentik ACM, Pearson SCP: The nature of volcanism. In Volcanic and Tectonic Hazard Assessment for Nuclear Facilities. Edited by: Connor CB, Chapman NA, Connor LJ. Cambridge University Press; 2009:74–115.View ArticleGoogle Scholar
- Costa A, Macedonio G: Computational modeling of lava flows: A review. Geological Society of America Special Papers 2005, 396: 209–218.Google Scholar
- Del Negro C, Fortuna L, Vicari A: Modelling lava flows by cellular nonlinear networks (CNN): preliminary results. Nonlinear Processes in Geophysics 2005, 12: 505–513. 10.5194/npg-12-505-2005View ArticleGoogle Scholar
- Diggle P: A kernel method for smoothing point process data. Applied Statistics 1985, 34: 138–147. 10.2307/2347366View ArticleGoogle Scholar
- Dragoni M: A dynamical model of lava flows cooling by radiation. Bulletin of Volcanology 1989, 51: 88–95. 10.1007/BF01081978View ArticleGoogle Scholar
- Dragoni M, Tallarico A: The effect of crystallization on the rheology and dynamics of lava flows. Journal of Volcanology and Geothermal Research 1994, 59: 241–252. 10.1016/0377-0273(94)90098-1View ArticleGoogle Scholar
- Duong T: Kernel density estimation and kernel discriminant analysis for multivariate data in R. Journal of Statistical Software 2007,21(7):1–16.Google Scholar
- Duong T, Hazelton ML: Plug-in bandwidth selectors for bivariate kernel density estimation. Journal of Nonparametric Statistics 2003, 15: 17–30. 10.1080/10485250306039View ArticleGoogle Scholar
- Griffiths RW: The dynamics of lava flows. Annual Review of Fluid Mechanics 2000, 32: 477–518. 10.1146/annurev.fluid.32.1.477View ArticleGoogle Scholar
- Harris AJL, Rowland SK: FLOWGO: a kinematic thermo-rheological model for lava flowing in a channel. Bulletin of Volcanology 2001, 63: 20–44. 10.1007/s004450000120View ArticleGoogle Scholar
- Harris AJL, Rowland SK: Controls on lava flow length. In Studies in Volcanology, The Legacy of George Walker, Special Publications for IAVCEI No.2, The Geological Society, London Edited by: Thordarson T, Self S, Larsen G, Rowland SK, Höskuldsson A. 2009, 33–51.Google Scholar
- Hill BE, Aspinall WP, Connor CB, Godoy AR, Komorowski JC, Nakada S: Recommendations for assessing volcanic hazards at sites of nuclear installations. In Volcanic and Tectonic Hazard Assessment for Nuclear Facilities. Edited by: Connor CB, Chapman NA, Connor LJ. Cambridge University Press; 2009:566–592.View ArticleGoogle Scholar
- Hornik K: The R FAQ.2009. [http://CRAN.R-project.org/doc/FAQ/R-FAQ.html] ISBN 3-900051-08-9Google Scholar
- IAEA: Volcanic hazards in site evaluation for nuclear power plants. In Draft Safety Guide DS405. International Atomic Energy Agency, Vienna, Austria; 2011.Google Scholar
- Jaquet O, Connor CB, Connor LJ: Probabilistic modeling for long-term assessment of volcanic hazards. Nuclear Technology 2008,163(1):180–189.Google Scholar
- Jarvis A, Reuter HI, Nelson A, Guevara E: Hole-filled seamless SRTM data. International Centre for Tropical Agriculture (CIAT) v4 edition. 2008. available from http://srtm.csi.cgiar.org available fromGoogle Scholar
- Kilburn CJR, Lopes RMC: The growth of aa lava fields on Mt. Etna, Sicily. Journal of Geophysical Research 1988, 93: 14,759–14,722.View ArticleGoogle Scholar
- Kilburn CJR, Luongo G: Active Lavas: Monitoring and Modelling. UCL Press, London, UK; 1993.Google Scholar
- Kiyosugi K, Connor CB, Zhao D, Connor LJ, Tanaka K: Relationships between temporal-spatial distribution of monogenetic volcanoes, crustal structure, and mantle velocity anomalies: An example from the Abu monogenetic volcano group, southwest Japan. Bulletin of Volcanology 2010., 71: doi:10.1007/s00,445–009–0316–4 doi:10.1007/s00,445-009-0316-4Google Scholar
- Le Bas M, Le Maitre RW, Woolley AR: The construction of the total Alkali-Silica chemical classification of the volcanic rocks. Mineralogy and Petrology 1986,46(1):1–22.View ArticleGoogle Scholar
- Luhr J, Simkin T: Paricutin: the Volcano Born in a Mexican Cornfield. Geocsience Press, Inc, Phoenix, AZ 1993.Google Scholar
- Malin MC: Lengths of Hawaiian lava flows. Geology 1980, 8: 306–308. 10.1130/0091-7613(1980)8<306:LOHLF>2.0.CO;2View ArticleGoogle Scholar
- Martin AJ, Umeda K, Connor CB, Weller JN, Zhao D, Takahashi M: Modeling long-term volcanic hazards through Bayesian inference: example from the Tohoku volcanic arc, Japan. Journal of Geophysical Research 2004.,109(B10208): doi:10.1029/2004JB003,201 doi:10.1029/2004JB003,201Google Scholar
- Miyamoto N, Sasaki S: Simulating lava flows by an improved cellular automata method. Computers and Geosciences 1997, 23: 283–292. 10.1016/S0098-3004(96)00089-1View ArticleGoogle Scholar
- Molina S, Lindholm CD, Bungum H: Probabilistic seismic hazard analysis: Zoning free versus zoning methodology. Bollettino di Geofisica Teorica et Applicata 2001,42(1–2):19–39.Google Scholar
- Pinkerton H, Wilson L: Factors controlling the lengths of channel-fed lava flows. Bulletin of Volcanology 1994, 56: 108–120.View ArticleGoogle Scholar
- Rowland SK, Gabriel H, Harris AJL: Lengths and hazards of channel-fed lava flows on Mauna Loa (Hawaii), determined from thermal and downslope modeling with FLOWGO. Bulletin of Volcanology 2005, 67: 634–647. 10.1007/s00445-004-0399-xView ArticleGoogle Scholar
- Silverman B: Density Estimation for Statistics and Data Analysis. No. 26 in Monographs on Statistics and Applied Probability, Chapman and Hall 1986.Google Scholar
- Silverman BW: Choosing the window width when estimating a density. Biometrika 1978, 65: 1–11. 10.1093/biomet/65.1.1View ArticleGoogle Scholar
- Stasiuk MV, Jaupart C: Lava flow shapes and dimensions as reflections of magma system conditions. Journal of Volcanology and Geothermal Research 1997,78(1–2):31–50. 10.1016/S0377-0273(97)00002-4View ArticleGoogle Scholar
- Stock C, Smith EGC: Comparison of seismicity models generated by different kernel estimations. Bulletin of the Seismological Society of America 2002, 92: 913–922. 10.1785/0120000069View ArticleGoogle Scholar
- Thordarson T, Self S: The Laki (Skaftar fires) and Grimsvotn eruptions in 1783–85. Bulletin of Volcanology 1993, 55: 233–263. 10.1007/BF00624353View ArticleGoogle Scholar
- Volentik ACM, Connor CB, Connor LJ, Bonadonna C: Aspects of volcanic hazard assessment for the Bataan nuclear power plant, Luzon Peninsula, Philippines. In Volcanic and Tectonic Hazard Assessment for Nuclear Facilities. Edited by: Connor CB, Chapman NA, Connor LJ. Cambridge University Press; 2009:229–256.View ArticleGoogle Scholar
- Wadge G, Young PAV, McKendrick IJ: Mapping lava flow hazards using computer simulation. Journal of Geophysical Research 1994, 99: 489–504. 10.1029/93JB01561View ArticleGoogle Scholar
- Walker GPL: Lengths of lava flows. Philosophical Transactions of the Royal Society, London 1973, 274: 107–118.View ArticleGoogle Scholar
- Wand MP, Jones MC: Kernel Smoothing. No. 60 in Monographs on Statistics and Applied Probability, Chapman and Hall 1995.Google Scholar
- Ward SN: A multidisciplinary approach to seismic hazard in southern California. Bulletin of the Seismological Society of America 1994, 84: 1293–1309.Google Scholar
- Woo G: Kernel estimation methods for seismic hazard area source modelling. Bulletin of the Seismological Society of America 1996, 88: 353–362.Google Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.