Preliminary assessment for the use of VORIS as a tool for rapid lava flow simulation at Goma Volcano Observatory , Democratic Republic of the Congo

Assessment and management of volcanic risk are important scientific, economic, and political issues, especially in densely populated areas threatened by volcanoes. The Virunga volcanic province in the Democratic Republic of the Congo, with over 1 million inhabitants, has to cope permanently with the threat posed by the active Nyamulagira and Nyiragongo volcanoes. During the past century, Nyamulagira erupted at intervals of 1–4 years – mostly in the form of lava flows – at least 30 times. Its summit and flank eruptions lasted for periods of a few days up to more than 2 years, and produced lava flows sometimes reaching distances of over 20 km from the volcano. Though most of the lava flows did not reach urban areas, only impacting the forests of the endangered Virunga National Park, some of them related to distal flank eruptions affected villages and roads. In order to identify a useful tool for lava flow hazard assessment at Goma Volcano Observatory (GVO), we tested VORIS 2.0.1 (Felpeto et al., 2007), a freely available software (http://www.gvb-csic.es) based on a probabilistic model that considers topography as the main parameter controlling the lava flow propagation. We tested different parameters and digital elevation models (DEM) – SRTM1, SRTM3, and ASTER GDEM – to evaluate the sensitivity of the models to changes in input parameters of VORIS 2.0.1. Simulations were tested against the known lava flows and topography from the 2010 Nyamulagira eruption. The results obtained show that VORIS 2.0.1 is a quick, easy-to-use tool for simulating lava-flow eruptions and replicates to a high degree of accuracy the eruptions tested when input parameters are appropriately chosen. In practice, these results will be used by GVO to calibrate VORIS for lava flow path forecasting during new eruptions, hence contributing to a better volcanic crisis management.


Introduction
During the past century, Nyamulagira (or Nyamuragira), the westernmost volcano of the Virunga Volcanic Province (VVP), erupted at intervals of 1-4 years (Burt et al., 1994;Smets et al., 2010aSmets et al., , 2015)).It is considered as one of the most active African volcanoes, with at least 39 documented eruptions since 1882 (Smets et al., 2010a).Most of historical eruptions of Nyamulagira occurred along its flanks, sometimes more than 10 km far from the central edifice.Nyamulagira eruptions commonly last few days to few weeks, but some voluminous and less frequent historical events lasted several months or years (Pouclet, 1976;Smets et al., A. M. Syavulisembo et al.: Rapid lava flow simulation at Goma Volcano Observatory 2015).Nyamulagira's activity is characterised essentially by lava fountaining activity along an eruptive fissure, which produces lava flows and progressively build, together with ejected tephra, a spatter-and-scoria cone along the fissure (Pouclet, 1976;Smets et al., 2015).
On 2 January 2010, an eruption of Nyamulagira started in the central caldera and along its SSE flank, emitting ∼ 45.5 × 10 6 m 3 of lava (Smets et al., 2014a).The related lava flows were of low viscosity, causing loss of vegetation in the Virunga National Park.The dense plume of gases that escaped from the eruptive vents affected the surrounding population (Cuoco et al., 2013).The estimated surface area covered by the 2010 lava flow is 15.17 ± 2.53 × 10 6 m 2 (Smets et al., 2014a).
Due to the presence of densely populated zones and a National Park in the vicinity of Nyamulagira, lava flow invasion probability associated with its frequent fissural eruptions must be assessed.In addition to detailed field surveys and literature reviews performed to better characterise Nyamulagira's eruptive activity and, hence, better assess eruption hazards, simulations are practical tools to create hazard maps and detect the most probably affected areas during an eruption.For lava flows, a wide diversity of both probabilistic and deterministic simulation models exist in the volcanological literature (e.g.Crisci et al., 1986Crisci et al., , 1997;;Ishihara et al., 1989;Wadge et al., 1994;Kuauhicaua et al., 1995;Felpeto et al., 2001;Favalli et al., 2005;Damiani et al., 2006) and offer different degrees of accuracy, depending on the mathematical methods used to make calculations and the input parameters considered.Unfortunately, most of these models only exist in the scientific literature and their source codes are not freely available.Furthermore, some require complex calculations and significant CPU usage, which is not always available.Felpeto et al. (2007) have built a volcanic hazard assessment software package (VORIS 2.0.1)within ArcGIS 9.1 ( © ESRI), which includes a probabilistic lava-flow model based on a previously developed model (Felpeto et al., 2001).Unlike existing lava-flow models, VORIS is free (downloadable from http://www.gvb-csic.es),easy to use and does not entail high computing requirements.However, VORIS may be less accurate than other more sophisticated deterministic lava-flow simulation models, which thus creates a dilemma as to whether a high precision or a quicker, easier-to-usebut less precise -tool is preferable.
The Goma Volcano Observatory (GVO) is in charge of conducting volcano monitoring and assessing volcanic hazards in the Democratic Republic of the Congo (DRC).Thus, as part of its work, it needs to adopt adequate methodologies and tools -tried and tested in other areas -that can be used by its scientists and technicians to comply with the most essential of its tasks and according to its scientific and technological resources.To check whether VORIS could be of use to GVO for lava flow hazard assessment in the DRC volcanoes, and for quickly estimating the potential impact of new lava flows in the event of a new eruption, we tested this package's lava-flow simulation model by replicating the Nyamulagira 2010 eruption, most of whose parameters and pre-eruption topography are known.
Nyamulagira (1.41 • S, 29.20 • E) is located in the Virunga National Park, which was declared as a World Heritage Site in 1979 and has an endangered one since 1994 (http://whc.unesco.org/fr/list/63).Nyamulagira lavas are mostly basic, SiO 2 -undersaturated and K-rich (Kampunzu et al., 1982;Aoki and Yoshida, 1983;Aoki et al., 1985), which form low viscosity flows able to travel for tens of kilometres.In recent decades, the eruptive activity took the form of periodic fissure-fed effusive eruptions, sometimes accompanied with summital activity (e.g.Hamaguchi and Zana, 1983;Wadge and Burt, 2011;Smets et al., 2014aSmets et al., , 2015)).Since April 2014, activity at Nyamulagira is restricted to episodic lava fountaining and the presence of a semi-permanent lava lake in a pit crater located in the NE sector of the central caldera (Smets et al., 2014b).

Methodology
We used VORIS 2.0.1 (Felpeto, 2009; http://www.gvb-csic.es) to simulate the lava flows emitted by Nyamulagira during the 2010 eruption.VORIS (Volcanic Risk Information Systems) was developed in an ArcGis 9.1 ( © ESRI) Geographical Information System (GIS) framework and allows eruptive scenarios and probabilistic hazards maps to be created rapidly for the commonest volcanic hazards (lava flows, fallout and pyroclastic density currents).
The lava flow model included in VORIS is a probabilistic (Monte Carlo algorithm) model based on the assumption that both topography and flow thickness play major roles in determining the path followed by a lava flow (Felpeto et al., 2007 and references therein).The model computes several possible paths for the flow on the basis of two simple rules: (i) the flow will only propagate from one cell to one of its eight neighbours if the difference in corrected topographic height between them is positive, and (ii) the probability that the flow will move from one cell to one of its neighbours is proportional to the difference in height.The calculation of the probability of a point being invaded by lava is performed by computing several random paths using a Monte Carlo algorithm (see Felpeto, 2009 for more details).
The input data for this simulation include a digital elevation model (DEM), from which we can chose either a single vent or a vent area including several vents, the maximum flow length (i.e. the total length of the path followed by the lava flow, not just the maximum longitudinal distance), and a height correction (i.e. the average thickness of the flow).All these input parameters remain essentially constant during the simulation and characterise the resulting model.The maximum flow length and the height correction or average thickness are integers defined in metres.In addition, the "iterations number" must be considered.This number is proportional to the number of calculations made and will determine the CPU time required and the accuracy of the results obtained.However, if a greater number of iterations suggests a more accurate results, it also increases the risk to overload the microprocessor.This is why in our study we paid special attention to this input parameter, in order to determine the most appropriate number of iterations for guaranteeing reliable results without using too much computing time.
VORIS allows to select one single vent (a single pixel) or a vent area (several pixels, e.g. a linear fissure).The Nyamulagira 2010 eruption occurred along an ENE-WSW-trending fracture in the S part of the caldera and a 600 m long fissure that opened on the SSE flank of the main edifice, which generated the lava flows we have simulated in this study.We simulated the vent area both as a single vent and as a fissure, obtaining very similar results due to the radial orientation of the fissure and the strong topographic control on lava emplacement.For simplification we only show results from a single vent, which we located at the middle of the eruptive fissure.
An accurate DEM is a fundamental requirement for modelling volcanic processes and associated risks, including pyroclastic, lava, and mud flows (Sheridan et al., 2004).In the present case study, the DEM plays a major role in determining the pathways of lava flows (Felpeto et al., 2007), as it is the numerical base used to compute the path the lava will follow, based on the principle used by VORIS.To test the model, we used three different DEMs, the best available ones to the GOMA Volcano Observatory scientists at the time of writing: SRTM1, SRTM3, and ASTER GDEM (Mukerejee et al., 2013).These three DEMs were created before the 2010 eruption, but only the ASTER DEM contains the topography of the neighbouring Nyamulagira 2006 eruption.The three DEMs also differ either by their resolution or by the way they were produced.SRTM (NASA's Shuttle Radar Topography Mission) DEMs were produced by radar interferometry, using radar images acquired during a space mission, with NASA's Space Shuttle.SRTM1 in the studied area has a spatial resolution of 31 m, while SRTM3 has 93 m-wide pixels.In recent years, SRTM DEMs were the main source of height data for the VVP.The ASTER GDEM (ASTER Global Digital Elevation Model) is obtained by stereophotogrammetry, using nadir and backward images acquired by the ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) optical sensor, onboard the TERRA satellite (NASA/METI/J-Spacesystems).* Fitness index to quantify the variation in the result of different simulations (calculated for all simulation using the fitness function of Favalli et al., 2009;Bilotta et al., 2012).
ASTER GDEM data are published with a spatial resolution of 30 m.
In order to evaluate the sensitivity of models to changes in input parameters, we first tested separately the influence of the flow length, the number of iterations and the lava-flow thicknesses using the SRTM1 DEM.We then conducted additional simulations using different DEMs with the most appropriate length, thickness and number of iterations determined by the first tests.The comparison of results obtained with the three different DEMs is important for testing the suitability of each model in this type of lava-flow simulations and for testing their aptness for use with VORIS.In equatorial zones, such as in the VVP, the quality of the ASTER GDEM is considerably affected by cloud cover.As a consequence, artifacts appear in the lava fields (Arefi and Reinartz, 2011).Compared to differential GPS ground control points acquired in the southern Nyiragongo lava field (Albino et al., 2015), the SRTM DEM has a better vertical accuracy than the ASTER GDEM.Hence, it is expected to have more reliable results using the SRTM DEMs.

Length parameter
The influence of the length parameter on the simulated flows was evaluated using SRTM1 and 12 different input length A vent location was fixed at point A (35N 746161/9840963) with an average lava flow thickness of 3 m.5000 iterations were applied with a different length value for each.Results are given in Figs. 2 and 3, and in with the values for the following parameters: 1. "True lava flow pixels": number of pixels of the DEM actually covered by the lava flow we are trying to simulate;    5,10,15,20,25,30,35,40,45,50, and 55 km length, respectively).
2. "Simulated pixels": number of pixels of the DEM corresponding to the surface covered by the simulated lava flow; 3. "Well-estimated pixels": number of pixels of the DEM corresponding to probable pixels that coincide with true pixels; 4. "Underestimated pixels": number of pixels of the DEM corresponding to true pixels that were not covered by probable pixels;  The results obtained reveal that an overly small length value underestimates the probability of being covered by the lava flow (simulation L1), while an overly large value (simulation L12) tends to overestimate the maximum longitudinal or run-out distance, as well as the lateral extent of the lava flow, given that the real eruption will stop before reaching this total maximum length.However, simulations L7 and L8, which considered intermediate maximum lengths of 30 and 35 km, respectively, match well the extent and run-out distance of the Nyamulagira 3 January 2010 lava flow.This indicates that, in this case and with the SRTM1 DEM, the total length of this lava flow was about double that of its longitudinal distance, thereby revealing the strong control that was exerted by the highly irregular topography characterised by non-rectilinear ravines and gullies.

Number of iterations
We arbitrary selected 10 different numbers of iterations (i.e.I1 = 10, I2 = 50, I3 = 100, I4 = 500, I5 = 1000, I6 = 5000, I7 = 10 000, I8 = 50 000, I9 = 100 000, I10 = 500 000).Simulations were carried out on SRTM1, with the same vent point A (UTM 35N, 746161 E/9840963 N), which corresponds to the flank vent location of the 2010 eruption, and fixed average thickness and total length of 3 m and 33 km, respectively.The results (Figs. 4 and 5 and   Increasing it further up to 8 m reduces the maximum distance reached by the lava flow but increases the surface covered by the simulated lava flow.This widening effect of the simulated flow is due to the fact that thicker flow can run over higher topography.As the path length is fixed at 33 km for all simulations, the runout distance of the simulated lava flow may consequently decrease for higher h if the topography is such that the lava flow may spread laterally, as it happens in the case simulated here.The result that best fits the 2010 lava flow is h3, which corresponds to a lava thickness of 3 m.that the number of simulated pixels, well-estimated pixels, and pixels outside the 2010 lava flow increase from simulations I1 to I6 but remain more or less constant from I6 to I10.Simulation I6 has the largest number of well-estimated pixels, the lowest number of underestimated pixels, and a relatively low running (calculation) time.Even so, it still gives over 15 000 pixels that are outside the actual lava flow.The number of iterations has less influence on the modelled lava flow length, with values ranging from 9.6 and 11.8 km.

Lava-flow thickness
In order to test the influence of the height correction parameter (average thickness) we performed eight simulations with different values of this parameter (i.e.h1=1 m, h2=2 m, h3=3 m, h4=4 m, h5=5 m, h6=6 m, h7=7 m, h8=8 m).Iterations and total lengths were fixed to 5000 and 33 km respectively.The results are presented in Figs. 6 and 7 and Table 3.
In these simulations, we notice that when the value of "h" increases, the effects of the relief decrease.In simulations h1-h8, the small zones enclosed by the true pixels decrease gradually in size until they are eliminated.Simulation h1 is too narrow and does not cover all the actual lava flow, while simulation h8 produces an overflow outside the 2010 lava flow.With the length parameter fixed to 33 km, we notice that the modelled lava-flow length increases for simulations h1 and h2.Simulation h3 has the greatest number of simulated and well-estimated pixels, but the lowest number of underestimated pixels.4. Simulation using SRTM1 (S1) gives the best match between the modelled and natural lava flows, as it corresponds to the best calibration determined in the previous sections.In simulation using SRTM3 (S2), lava  spreads over a greater surface area and clearly exceeds the real area covered by the real lava flow.This result highlights the need to calibrate input parameters according to the DEM used in the simulation.The simulation S3 (ASTER GDEM) stopped prematurely because of strong artifacts present in the DEM.ASTER GDEM seems consequently not appropriate for lava flow modelling in the VVP.
It is worth mentioning that we were aware of the existence of other DEMs with higher resolutions.However, computational requirements to download and process them with VORIS are not available at GVO at the time of writing.

Discussion and conclusions
We conducted several simulations of lava-flow emplacement using VORIS 2.0.1 (Felpeto et al., 2007), in order to replicate the lava flow emitted by the 2010 Nyamulagira eruption and test the aptness of this model for conducting lava-flow hazard assessment at the Nyamulagira volcano.Model calibration was performed by changing the values of the model parameters until simulations best matched the 2010 lava flow.Next, this calibrated simulation was performed on three different DEMs to detect their influence on results.
This analysis shows that some input parameters can drastically change results.Too few iterations produce very poor results with a high degree of inaccuracy, while too many lead to considerable overestimates and a consequent increase in computing time.In our study, a number of iterations between 5000 and 10 000 gave the best results and the greatest de- gree of coincidence (see Tables 1-4 and Figs . 3, 5, 7, and 9) between the simulations and the real lava flow, with a total computing time of less than 20 min for each run.Simulations were less sensitive to changes in height correction (i.e.average thickness).However, we observed that the best results were obtained with a value of 3 m, which coincides with the average thickness of the actual lava flow.Concerning the length parameter, our results indicate that an overly small length value underestimates the probability of being covered by the lava flow, while an overly large value tends to overestimate the maximum longitudinal or run-out distance, as well as the lateral extent of the lava flow.This highlights the strong influence of irregular topography on the final coverage estimates.Results also show that the model calibration must be adapted to the DEM used, as both spatial resolution and DEM quality strongly affect the results.
In summary, models generated using VORIS 2.0.1 succeeded to replicate the emplacement of the Nyamulagira 2010 lava flow with a reasonable degree of accuracy.Considering that this is a quick, easy-to-use software and one of the few that are freely available on the Internet, it is perfectly adapted for lava flow hazard assessment performed at Goma Volcano Observatory.Further work will be dedicated to establish a method for determining the best input parameters (lava flow length, thickness and number of iterations) to be used when no a priori lava flow surface is known.This will be done by testing the modelling with other known lava flows from Nyamulagira.Additional work will also be required to adapt simulations to the neighbouring Nyiragongo volcano, for which the flank eruption style and lava flow propagation are different.

Figure 1 .
Figure 1.Location of the Virunga region with the Nyamulagira volcano and 2010 lava flows indicated in red.

Figure 2 .
Figure 2. Probability of lava flows invasion from the emission point A for various lengths values (L1-12).As for the rest of Figs showing the results of simulations, we show here the simulated areas corresponding to a hazard > 0.1 % when simulated result matches the real lava (well-estimated pixels), and to a hazard < 0.1 % when the simulated result falls outside the real lava flow (outside estimated pixels).

Figure 6 .
Figure 6.Simulations considering different values of height correction ( h) (average lava flow thickness).Increasing the lava flow thickness from 1 to 3 m increases the maximum distance reached by the lava flow.Increasing it further up to 8 m reduces the maximum distance reached by the lava flow but increases the surface covered by the simulated lava flow.This widening effect of the simulated flow is due to the fact that thicker flow can run over higher topography.As the path length is fixed at 33 km for all simulations, the runout distance of the simulated lava flow may consequently decrease for higher h if the topography is such that the lava flow may spread laterally, as it happens in the case simulated here.The result that best fits the 2010 lava flow is h3, which corresponds to a lava thickness of 3 m.
conducted several simulations to test the suitability of the different DEMs (S1: SRTM1, S2: SRTM3, and S3: ASTER GDEM) using vent location at point A (UTM 35N, 746161 E/9840963 N), lava thickness of 3 m, 5000 iterations, and 33 km for the length parameter.Results are shown in Figs. 8 and 9, and in Table

Figure 8 .
Figure 8. Probability of lava flows invasion from A (emission point of the 2010 eruption of Nyamulagira) for different DEMs (S1, S2, and S3).

Table 1 .
Influence of the length parameter on the simulation of the 2010 lava flow of Nyamulagira.

Table 2 .
Influence of the number of iterations used for the simulation of the 2010 lava flow of Nyamulagira.

Table 3 .
Influence of lava flow average thickness on the simulation of the 2010 lava flow of Nyamulagira.

Table 4 .
Influence of the DEM used for the simulation of the 2010 lava flow of Nyamulagira.