Assessing lava flow susceptibility at neighbouring volcanoes: Nyamulagira and Nyiragongo volcanoes, Virunga Volcanic Province

Abstract

Assessing volcanic hazards in locations exposed to multiple central volcanoes requires to consider multiple potential eruption sources and their respective characteristics. While this is common practice in ashfall hazard assessment, this is generally not considered for topography-controlled volcanic flow processes. Yet, in volcanic areas with closely spaced volcanic systems, eruptions fed from several contrasted volcanic systems might threaten one given area. Considering the case of the Nyiragongo and Nyamulagira volcanoes in the Virunga Volcanic Province (D.R.Congo), we present a method to produce a combined lava flow inundation susceptibility map that integrates both volcanoes. The spatial distribution of the probability of vent opening for the next eruption is separately constrained for both volcanoes based on the mapping of historical and pre-historical eruptive vents and fissures. The Q-LavHa lava flow probability model is then calibrated separately for each volcano, considering several historical lava flows of Nyamulagira (2004, 2006, 2010) and Nyiragongo (2002). The maps for the two volcanoes are thereafter integrated based on a weighted sum of both individual lava flow inundation probability maps, assuming historically-based relative eruption frequency of the two volcanoes. The accuracy of this probabilistic susceptibility map for the most active volcanic region in Africa was unfortunately validated by the May 2021 lava flow produced by Nyiragongo. This map was discussed and validated in 2019 with local scientists, as well as representatives of disaster management and urban planning institutions, but was not included in the regional contingency plan ahead of the 2021 eruption crisis. Updating the volcanic crisis and evacuation management plans with this lava flow probability map could contribute to reinforce risk awareness among the population and inform the future development of the city of Goma.

Introduction

Assessing volcanic hazards and risks in locations exposed to multiple volcanic sources is challenging. It requires combining the multiple hazards and their characteristics (frequency, magnitude) from each source individually to assess the cumulative hazard (and associated risk) at any given location. This challenge has been already tackled in highly active volcanic regions subjected to, among others, ash fallout and pyroclastic flows from multiple stratovolcanoes (e.g. Lirer et al. 2010; Hurst and Smith 2010; Alberico et al. 2011; Jenkins et al. 2012; Biass et al. 2014; Bevilacqua et al. 2015; Neri et al. 2015; Jiménez et al. 2020).

In the context of long-term lava flow hazard assessment, multiple sources of eruptions are generally considered within a monogenetic volcanic field or on a polygenetic volcano with scattered parasitic vents (e.g. Gallant et al. 2018). These studies commonly approach the spatial probability of future vent opening through the estimation of the eruption susceptibility (e.g. Bartolini et al. 2013; Becerril et al. 2013; Cappello et al. 2012; Marti and Felpeto 2010). Lava flow models are subsequently used to constrain the probability of lava flow inundation, assuming similar or gradually varying lava flow characteristics (e.g. lava flow length, geometric properties, viscosity and effusion rate) from any location with a positive eruption susceptibility (e.g. Favalli et al. 2011; Favalli et al. 2012; Cappello et al. 2015; Del Negro et al. 2019). Assessing lava flow hazards is more challenging when one geographical area is exposed to potential inundation by lava flows from contrasted volcanic systems characterised by very different lava flow characteristics and eruption frequency. This is the case in the proximity of closely spaced effusive volcanoes, e.g. Mauna Loa – Kilauea (United States), Askja – Bárðarbunga – Grimsvötn (Iceland), Katla – Eyjafjallajökull (Iceland), or Nyamulagira – Nyiragongo (Democratic Republic of Congo).

Nyamulagira and Nyiragongo are a particular case in this regard since these very active volcanoes are located less than 15 km from each other and are characterised by highly contrasted lava flow characteristics, including viscosity, lava flow length and eruption frequency (Morrison et al. 2020; Smets et al. 2015). Nyamulagira and Nyiragongo also represent important risk implications, as both volcanoes are located in the close proximity of the densely populated area along the northern shoreline of Lake Kivu, including the more than one million inhabitants of the urban areas of Goma (D.R.C.) and Gisenyi (Rwanda). These areas were heavily affected by lava flows and ground fracturing during the 2002 and 2021 eruptions of Nyiragongo volcano (Fig. 1, Allard et al. 2002; Komorowski et al. 2002; Michellier et al. 2020).

Fig. 1

Oblique 3D view looking North of the city of Goma and Nyiragongo and Nyamulagira volcanoes produced using Google Earth. The foreground shows the city of Goma bordering Lake Kivu, with highlights of Mt Goma volcanic cone (light green) and the paths of the three historical lava flow eruptions. Only the historical lava flows of Nyiragongo volcano are highlighted in red colours. Goma airport and the border with Rwanda are visible to the East

Previous studies assessed the lava flow hazard at Nyiragongo volcano (Chirico et al. 2009; Favalli et al. 2009; Copernicus 2018) and modelled lava flow paths through the city of Goma (Favalli et al. 2006). These studies demonstrated that the areas destroyed by the 2002 and 2021 lava flows are characterised by a high susceptibility to lava flow inundation, whereas other areas are characterised by a much lower susceptibility. The focus then was largely on the modelling of lava flows originating from the southern flank of Nyiragongo volcano and the potential impact of these lava flows on the city of Goma. The lava flow hazard originating from other parts of the volcanic field and the consideration that some areas could potentially be exposed to lava flow from both Nyamulagira and Nyiragongo has not yet been investigated.

We present a new integrated lava flow susceptibility map for Nyamulagira and Nyiragongo volcanoes. We refer to ‘susceptibility map’ and not to ‘hazard map’ as the maps present the relative spatial variation in the probability of lava flow inundation, without constraining this probability within a given time period. Individual lava flow susceptibility maps for Nyamulagira and Nyiragongo volcanoes are calibrated, accounting for the contrasted vent opening susceptibility and lava flow characteristics (i.e. length and thickness) of these two neighbouring volcanoes. We investigate two potential scenarios to combine the lava flow susceptibility of Nyamulagira and Nyiragongo to produce a combined lava flow susceptibility map for both volcanoes. We then consider the implication of this new map for the assessment of risk for the densely populated urban areas along the northern shoreline of Lake Kivu and explain how this map was presented and received by provincial authorities of North Kivu, prior to the 2021 eruption crisis.

Volcanic activity in the Virunga Volcanic Province

The Virunga Volcanic Province (VVP) is located within the western branch of the East African Rift System, at the intersection of the Democratic Republic of Congo (D.R.C.), Rwanda and Uganda (Fig. 2). The VVP comprises eight large volcanic edifices, including the active Nyamulagira and Nyiragongo volcanoes, and hundreds of pyroclastic cones scattered over the lava plain and flank of the Virunga volcanoes (Denaeyer 1975; Capaccioni et al. 2002; Poppe et al. 2016a). Historic volcanic activity in the VVP has been restricted to the western intra-rift part of the VVP, with the exception of the 1957 eruption of Mugogo volcano (Fig. 2; (Verhaeghe 1958).

Fig. 2

Location of Nyamulagira and Nyiragongo volcanoes. A Map of the Nyamulagira and Nyiragongo volcanic field showing the location of eruptive fissures and vents, and the historical lava flows of Nyamulagira and Nyiragongo volcanoes (from Smets et al. 2010; Smets and Poppe 2016). B Map showing the main volcanic edifices of the VVP (white triangles), the location of the southern part of the Virunga National Park (green area), roads (yellow lines) and urban areas (red areas). Nya = Nyamulagira; Nyi = Nyiragongo; Mik = Mikeno; Kar = Karisimbi; Vis = Visoke; Sab = Sabinyo; Gah = Gahinga; Muh = Muhavura. C Location of the VVP within the western branch of the East African Rift System. Background: SRTM DEM based hillshade image

Nyamulagira is one of the most active volcanoes on the African continent, with at least 42 documented lava flow eruptions since 1882 (one eruption every 1 to 4 years; (Smets et al. 2015). Over the period 1882 until 2014, the eruptive activity of Nyamulagira volcano was mainly characterized by lava flows issued from eruptive fissures, mostly along the NNW – SSE oriented fissure network which crosses the main edifice (Smets et al. 2015). The composition of lava flows from Nyamulagira volcano are alkali basalts, hawaiites, basanites and tephrites with a range of silica content from 43–56 wt.% (Morrison et al. 2020). Nyamulagira lava flows have higher viscosity and yield strength than Nyiragongo flows (Morrison et al. 2020): they are typically emplaced over several days to weeks, and have an average thickness of 3 m (Smets et al. 2015). Since 2014, the eruptive activity of Nyamulagira volcano is characterised by the intermittent presence of a lava lake in the pit crater of the central caldera (Coppola et al. 2016). The lava flows of Nyamulagira volcano represent a major threat to the Virunga National Park, an endangered UNESCO World Heritage site, and to the inhabitated area of Sake, that were affected four times by lava flows during the first half of the 20^th century (Fig. 2; Smets et al. 2014).

Nyiragongo volcano is renowned for the presence of a semi-permanent summit lava lake and the eruption of highly alkaline and strongly silica-undersaturated lavas (Platz 2002; Chakrabarti et al. 2009). The volcano produces highly fluid lava flows that move rapidly at flow velocity of tens of kilometres per hour (Tazieff 1977; Komorowski et al. 2002; Syavulisembo 2019). Over the last 120 years, the volcanic activity of Nyiragongo volcano has been characterised by intra-crater activity linked to the presence of the semi-permanent lava lake showing large level fluctuations (Durieux 2002a; Komorowski et al. 2002; Smets et al. 2016). This intra-crater activity was interrupted by flank eruptions in 1977, 2002 and 2021 (Fig. 1). During each of these events, a set of parallel fissures opened on the flanks of Nyiragongo volcano, feeding several lava flow branches of volumes of 10-25 ×10⁶ m³ in a few hours (Tazieff 1977; Cocheme and Vellutini 1979; Allard et al. 2002; Komorowski et al. 2002; Tedesco et al. 2007; Syavulisembo 2019; Smittarello et al. 2022). The rapidly advancing lava flows from each of these eruptions reached, to various extents, inhabited areas of Nyiragongo southern flank, destroying buildings and causing several tens to hundreds of fatalities (Tazieff 1977; Durieux 2002b; Smittarello et al. 2022). Further extension of the eruptive fissures to the south during the 2002 event caused two lava flow branches to cross the city of Goma, one of them reaching Lake Kivu, destroying about 10% of the urban area (Syavulisembo 2019). Both the 2002 and 2021 eruption crises caused the spontaneous evacuation of hundreds of thousands of Goma inhabitants ( Komorowski et al. 2002; Smittarello et al. 2022).

Methodology

The distinct characteristics of Nyamulagira and Nyiragongo volcanoes were accounted for by producing lava flow susceptibility maps separately for both volcanoes before combining them. The following sections detail the main methodological steps. First the vent opening susceptibility was modelled for each volcano using QVAST (Bartolini et al. 2013) based on mapped eruptive fissures and vents. An analysis of the length of historical lava flows was also performed. Based on these input data, lava flow model calibration and simulation of the lava flow susceptibility was conducted separately for each volcano using Q-LavHA (Mossoux et al. 2016). Both lava flow susceptibility maps were subsequently combined to produce a lava flow susceptibility map for the intra-rift part of the VVP (Fig. 3).

Fig. 3

Flowchart presenting the main steps of the methodology used to derive a combined lava flow susceptibility map for Nyiragongo and Nyamulagira volcanoes

Vent opening susceptibility

The vent opening susceptibilities at Nyamulagira and Nyiragongo volcanoes were calculated using the QGIS for Volcanic Susceptibility (QVAST) plugin developed by Bartolini et al. (2013). Vent opening susceptibility assumes that the spatial distribution of the probability of vent opening of the next eruption can be modelled by integrating the spatial distribution of different volcanic or magmatic features (e.g. vent, fissures, dykes; Cappello et al. 2012; Bartolini et al. 2013; Becerril et al. 2013). Here, the locations of eruptive fissures and vents were obtained from the volcano-structural map of Smets and Poppe (2016) and were attributed to the Nyamulagira or Nyiragongo volcanic fields based on the mapping of Thonnard et al. (1965). These features were then used as input volcano structural data in QVAST. No distinction was made between historical and non-historical fissures and vents in the development of the lava flow susceptibility map because the eruptive history of both volcanoes is only constrained since the start of the 20^th century and few historical vents can be characterized at Nyiragongo. For both volcanoes separately, probability density function (PDF) maps were derived by calculating a gaussian kernel density function on the eruptive fissures and vents. The Least Square Cross Validation (LSCV) criterion provided in QVAST was used to estimate the optimal bandwidth (h) of each kernel density function. The resulting PDF maps for fissures and vents of each volcano were combined through an equal weighted sum to provide the final vent opening susceptibility map for each volcano.

Lava flow model calibration

The Q-LavHA lava flow emplacement model (Mossoux et al. 2016) was used to constrain the lava flow inundation susceptibility. Q-LavHA is an open source probabilistic lava flow inundation model running as a plugin in QGIS and simulates lava flows sourced from a point, line or regularly spaced grid according to the probabilistic steepest slope (Mossoux et al. 2016). The model accounts for (1) the ability of lava flows to overcome small topographic obstacles and fill depressions through the use of two corrective factors (i.e., Hc and Hp, respectively) and (2) the vent opening susceptibility when simulating lava flows from a regularly spaced grid.

The optimal values of Hc and Hp corrective factors, which can be related respectively to the mean and maximum lava flow thickness, were determined for each volcano by simulating the extent of the 2004, 2006 and 2010 lava flows for Nyamulagira volcano, and the different lava flows of the 2002 eruption for Nyiragongo. These flows were selected because digital topographic data pre-dating their emplacement is available and because they represent a range of emplacement conditions (i.e. contrasted flank slope, eruption rates). As the lava flows considered for model calibration should not be present on the DEM, the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) acquired in 2000 at a 30-meter resolution was used as topographic base for model calibration. The initial Hc and Hp corrective factors of each lava flow branch of the 2002 eruption of Nyiragongo were based on documented lava flow thickness values (Favalli et al. 2006; Favalli et al. 2009). For the 2004, 2006 and 2010 Nyamulagira lava flows, no specific thickness values were available. The initial Hc corrective factor for these lava flows was set to 3 m, which corresponds to the average thickness of lava flows erupted by Nyamulagira volcano (Pouclet 1976; Smets et al. 2015). The Hc and Hp corrective factors were subsequently adapted through an iterative process until the simulated lava flow extent provided an optimal fit with the real lava flow obtained from the mapping of Smets and Poppe (2016). The latter was determined based on the true positive, false positive and false negative fitness indices provided in Q-LavHA (Mossoux et al. 2016), derived from comparing the agreement between the area of the simulated and the reference lava flow.

Length of historical lava flows

The length of the historical lava flow branches of Nyamulagira and Nyiragongo volcanoes was systematically measured on the volcano-structural map of Smets and Poppe (2016). For Nyiragongo volcano, all the lava flow branches of the 1977 and 2002 eruption were considered except the few short lava flows that occurred on the upper flank of Nyiragongo volcano above an elevation of 2,700 m because of their very small length compared to classical lower flank Nyiragongo lava flows that can threaten inhabited areas. For Nyamulagira volcano, the lava flows of the eruptions of 1938 – 1940, 1948, 1951 – 1952, 1954, 1956, 1957, 1958, 1967, 1971, 1976 – 1977, 1980, 1981 – 1982, 1984, 1986, 1987 – 1988, 1989, 1991 – 1993, 1994, 1996, 1998, 2000, 2001, 2002, 2004, 2006, 2010 and 2011 – 2012 were considered (Smets et al. 2010; Smets and Poppe 2016). Model calibration by Mossoux et al. (2016) highlighted that the lava flow length should be increased proportionally to the DEM resolution to simulate the entire extent of the lava flow, due to the meandering effect of individual iteration paths. A correction factor was derived from the model calibration of the 2002 lava flow branches of Nyiragongo volcano and the 2004, 2006 and 2010 lava flow branches of Nyamulagira volcano to correct for the impact of the DEM resolution on the simulated lava flow length. Afterwards, the distribution of the lava flow lengths for each volcano, multiplied by the obtained correction factor, was approximated by a normal distribution, truncated at 0 km length, considering the mean and standard deviation of measured lava flow lengths. The probability for the lava flow to reach a certain length was computed as the exceedance probability function derived from the normal distribution of lava flow length (Bonne et al. 2008).

Lava flow inundation susceptibility

The lava flow inundation susceptibility was calculated separately for Nyamulagira and Nyiragongo volcanoes in Q-LavHA by using the derived vent opening susceptibility maps, volcano-average Hc and Hp corrective factors and lava flow length statistical distribution. A TanDEM-X derived DEM (Albino et al. 2015), degraded to 30 m resolution to avoid local topography artefacts, was used as topographic base. This DEM was preferred to the SRTM, as it represent the topography after the last 2011-2012 lava flow eruption of Nyamulagira. Simulations were run before the occurrence of the 2021 Nyiragongo eruption. For a limited area on the western rift shoulder where no TanDEM-X derived DEM was available, the SRTM 30 m DEM was used as topographic base (Fig. 2).

One hundred lava flow paths were simulated from every pixel of the DEM (i.e. every 30 m) and weighted according to the vent opening susceptibility map. The value of 100 lava flow path iterations was derived from sensitivity tests assessing the influence of the vent spacing and number of iterations on model outputs. These tests showed that when lava flow paths are simulated from each pixel, the output of a lava flow simulation with one hundred or five hundred lava flow paths is equivalent, while hundred iterations significantly reduced the model processing time (Figs. 10, 11 and 12 in Appendix 1). This is so because many iterations from adjacent pixels channelize in the same topographic lows, leading to stability of model outputs. In total ca. 1.8×10⁸ and 1.0×10⁸ lava flow paths were simulated for Nyamulagira and Nyiragongo volcanoes, respectively.

The lava flow inundation susceptibility maps of Nyamulagira and Nyiragongo volcanoes express the susceptibility of each map pixel to lava flow inundation during the next eruption of each of those volcanoes (i.e. the sum of all the probabilities equals one). The probability prob_(i,j) of a pixel i to be inundated by lava flow from a vent j, is calculated as the sum of all the k path_(i,j) that effectively reach point i from vent j, weighted by the probability for the lava to reach the length that was simulated between the vent j and the pixel i following path k (${\overline{\mathrm{\varphi }} }_{{({\text{i}},{\text{j}})}_{k}}$), divided by the total number of paths n simulated from vent j (Eq. 1):

$${{prob}_{(i,j)}}=\frac{\sum_{k=0}^{{\text{n}}}({{\text{path}}}_{{({\text{i}},{\text{j}})}_{k}}\times {\overline{\mathrm{\varphi }} }_{{({\text{i}},{\text{j}})}_{k}})}{{\text{n}}}$$

(1)

Where path_(i,j) is a binary variable that is equal to one if the flow iteration from vent j effectively flow to pixel i, and is zero otherwise.

The susceptibility of lava flow inundation at each pixel i (S_i) is calculated as the sum for all potential vents j of the probability to inundate pixel i (prob_(i,j), equation1), weighted by the probability of vent opening at the vent j (p_ventj) divided by the total number of vents j_max simulated within the volcanic field (Eq. 2).

$${{\text{S}}}_{i}=\frac{\sum_{{\text{j}}=1}^{{j}_{max}}{({\text{prob}}}_{({\text{i}},{\text{j}})}{\times {\text{p}}}_{{vent}_{j}})}{{j}_{max}}$$

(2)

The values of S_i for all pixels vary between 0 and 1 within a lava flow field, but most values are very small (10^-17 to 10^-5) and non-normally distributed with a strong skewness towards higher values. As the sum of all susceptibility values in a volcanic field is equal to one, the valuse of individual pixel are also dependent on the size of the volcanic field. For the purpose of visual contrast and comparison of different susceptibility maps, they are displayed using percentile classes of susceptibility values on the maps. For instance, the pixels within the 98-100 percentile class correspond to the 2% of the pixels with the highest S_i value over the entire flow field. Up to this point, the lava flow susceptibility map of each volcano is independent from the timing of the next eruption and the relative frequency of eruptions at the two volcanoes.

Combination of the lava flow susceptibility maps

In order to produce a long-term lava flow susceptibility map for the intra-rift part of the VVP (S_VVP), the individual lava flow susceptibility maps of Nyamulagira (S_Nya) and Nyiragongo (S_Nyi) volcanoes were combined at pixel scale as follows (Eq. 3):

$${{\text{S}}}_{{\text{VVP}}}={{\text{W}}}_{{\text{Nya}}} *{{\text{S}}}_{{\text{Nya}}}+ {{\text{W}}}_{{\text{Nyi}}} * {{\text{S}}}_{{\text{Nyi}}}$$

(3)

where W_Nya and W_Nyi are weights that are related to frequency of lava flow eruptions at Nyamulagira and Nyiragongo volcanoes, with (W_Nya + W_Nyi = 1). As the long-term temporal probability of lava flow eruptions at Nyiragongo and Nyamulagira is poorly constrained, a reference scenario is produced in which the two susceptibility maps are combined using equal weight, assuming a similar temporal probability of eruption at Nyamulagira and Nyiragongo (i.e. W_Nya = W_Nyi = 0.5). This scenario is considered as reference, not because it is considered as most realistic, but because it makes no assumption on the relative frequency of eruption between the two volcanoes. Using the limited historical records since the mid-20^th century (i.e. the period for which the observations are fully exhaustive), a total of 30 independent lava flow eruptions have been documented at Nyamulagira (between 1948 and 2012, Smets et al. 2015) with only two eruptions being documented at Nyiragongo over the same period at the time of the simulation. Although this historical record is very short and does not allow derivation of an absolute eruption frequency at Nyiragongo, we develop a second scenario to integrate the two susceptibility models assuming that eruptions at Nyamulagira are 15 times more likely than those at Nyiragongo, based on the relative number of eruptions over 72 years (i.e. W_Nya = 0.937 and W_Nyi = 0.063). This scenario serves as a numerical experiment to illustrate the importance of considering eruption frequency of the different volcanic systems when integrating susceptibility maps from different volcanoes.

In order to analyse the difference in lava flow susceptibility between the combined and individual lava flow susceptibility models at any given location, the difference in lava flow susceptibility percentile value (${\Delta }_{percentile}$) was calculated for both scenarios as follows (Eq. 4):

$${\Delta }_{percentile}={SP}_{VVP}-{\text{max}}\left({SP}_{Nya},{SP}_{Nyi}\right)$$

(4)

Where ${SP}_{VVP}$ is the lava flow susceptibility percentile map for the intra-rift part of the VVP, and ${SP}_{Nya}$ and ${SP}_{Nyi}$ are the lava flow susceptibility percentile map of Nyamulagira and Nyiragongo volcanoes, respectively.

Results

Vent opening susceptibility

Figure 4 illustrates the long-term vent opening susceptibility for the next eruption of Nyamulagira and Nyiragongo volcanoes. The optimal bandwidths of the PDF maps of the fissures and vents at both volcanoes are approximately similar, and the bandwidths of the PDF maps of the fissures are about half of those of the vents for both volcanoes (Nyamulagira fissures: 1,056 m; Nyiragongo fissures: 1,230 m; Nyamulagira vents: 2,040 m; Nyiragongo vents: 2,106 m). The vent opening susceptibility at Nyamulagira volcano is the highest along the NNW – SSE oriented fissure system going across the volcano summit, and at the location of some distant flank eruptions (Fig. 4A) consistent with the patterns described by Smets et al. (2015). The vent opening susceptibility at Nyiragongo volcano is the highest along the N – S oriented fissure system on the southern flank that hosted the main 1977 and 2002 fissures, but also on the NW flank. The 2021 eruptive vents and fissures, not taken into account in producing this vent opening susceptibility map, opened in this highest-susceptibility zone on the upper to mid southern flank, but a minor fissure also opened at the NW base of Nyiragongo. Another significant alignment of vents and fissures – the Rushayo alignment described by Denaeyer (1975) - extends from the SW flank of Nyiragongo to the Lake Kivu shoreline, defining a secondary maximum of vent opening susceptibility in the southern western part of the Nyiragongo volcanic field (Fig. 4B). The vent opening susceptibility of both volcanoes shows an overlap of ca. 360 km2 at the border between the Nyamulagira and Nyiragongo volcanic fields (Fig. 4C). The vent opening susceptibility for both volcanoes is generally low to moderate in this overlapping zone, except for the southern Rushayo alignment.

Fig. 4

Long-term vent opening susceptibility during the next eruption of Nyamulagira (A) and Nyiragongo volcano (B). Map frame C shows the overlapping area between the vent opening susceptibility maps of Nyamulagira and Nyiragongo volcanoes. The dashed line shows the location where the vent opening susceptibility is equal to zero. Background: hillshade relief map of the SRTM 30 m DEM

Lava flow parameter calibration

Figure 5 shows the simulations for the 2004, 2006 and 2010 lava flows of Nyamulagira, and the lava flow branches of the 2002 eruption of Nyiragongo obtained with Q-LavHA. Table 1 provides the optimal parameters that were used to simulate these lava flows. For the 2004, 2006 and 2010 lava flows of Nyamulagira volcano, the Hc corrective factors range between 1.5 and 2 m, and the Hp corrective factors range between 2.5 and 4 m. These corrective factors are of the same order of magnitude as the average lava flow thickness observed at Nyamulagira volcano (i.e. 3 m; Smets et al. 2015). For the 2002 lava flow of Nyiragongo volcano, the Hc corrective factors range between 0.5 and 2 m, and the Hp corrective factors range between 1 and 3 m. These corrective factors are consistent with the lower viscosity and thickness of lava reported for Nyiragongo volcano (Morrison et al. 2020). Table 1 provides an evaluation of each simulated lava flow with the true positive being the proportion of the simulated lava flow field that corresponds to the actual lava flow outline, and the false negative, being the proportion of the actual lava flow with zero inundation susceptibility (Mossoux et al. 2016). In general, the lava flow emplacement model simulates well the main branches of the lava flows of Nyamulagira and Nyiragongo volcanoes, including their contrasted width, although true positive values remain moderate (Mossoux et al. 2016). A low true positive value and a high false positive value is obtained for the simulation of the 2006 lava flow of Nyamulagira because the simulated southernmost branch is much longer than the real branch. It is to be noted that all pixels with a non-zero probability are considered for the calculation of these accuracy parameters: in practice, high probability values are mainly constrained to the actual flow extent and false positives are mostly associated with low probability pixels. Based on the calibration of Q-LavHA for the different flows of each volcano, the Hc and Hp parameters were respectively set to 2 and 3 m for Nyamulagira volcano and 1 and 3 m for Nyiragongo volcano, the lower Hc of Nyiragongo generally inducing narrower flows.

Fig. 5

Lava flow simulations obtained with Q-LavHA for the different branches of the 2002 lava flow of Nyiragongo volcano (A – F) and the 2004, 2006 and 2010 lava flows of Nyamulagira volcano. See simulation parameters in Table 1

Table 1 Parameters of the lava flow simulations obtained with Q-LavHA for the 2004, 2006 and 2010 lava flows of Nyamulagira volcano, and the 2002 lava flow of Nyiragongo volcano. The names of the 2002 lava flow branches of Nyiragongo volcano refer to Fig. 5

Length of historical lava flows

Figure 6 shows the frequency of occurrence of lava flow lengths for Nyamulagira and Nyiragongo volcanoes and the decreasing cumulative frequency, which can be approximated with a normal distribution. These frequency distributions highlight the longer lava flow length of Nyamulagira compared to Nyiragongo volcano. The average and standard deviation of the length of the historical lava flow branches was estimated at 12,860 ± 5,980 m for Nyamulagira volcano, and 7,530 ± 2,050 m for Nyiragongo volcano. A similar average length correction factor of 1.66 was estimated for both Nyamulagira and Nyiragongo volcanoes from the sets of calibrated lava flows (Table 1) and was implemented in the simulation.

Fig. 6

Frequency of occurrence (grey bars), probability density distribution (grey line) and exceedance probability (dashed line) of lava flow length of Nyamulagira and Nyiragongo volcanoes. The very small lava flows that occurred on the upper flank of Nyiragongo volcano above an elevation of 2700 m were not included in the analysis of the lava flow length and are therefore not included in the figure

Long-term lava flow susceptibility

Figures 7 and 8 show the susceptibility of lava flow inundation during the next eruption for Nyamulagira and Nyiragongo volcanoes, respectively, assuming the vent opening susceptibility shown in Fig. 4. At Nyamulagira, high to very high probabilities of lava flow invasion are found on the lower NW and SE flanks, particularly in the areas that were already invaded by lava flows over the period 1948 – 2012 (Fig. 7B). The districts located in the western part of Goma and the city of Sake are characterized by an intermediate to low probability of lava flow invasion by Nyamulagira flows, but several segments of the Goma-Sake road have a high susceptibility to be invaded by lava. The villages located along the North-South road Goma-Rutshuru and in the proximity of the western rift shoulder are characterised by a low to intermediate probability of lava flow invasion. High to very high inundation susceptibilities are mainly constrained within the boundaries of the Virunga National Park, to the exception of the village of Kahe and its surroundings, north-west of the volcano (Fig. 7).

Fig. 7

Frame A illustrates the lava flow inundation susceptibility map for Nyamulagira volcano. B The areas that have been invaded by lava flows over the period 1948 – 2012 are superimposed on the lava flow susceptibility map in frame B. Class boundaries are defined by susceptibility percentiles. Absolute susceptibility values for each class boundary are as follows: 0^th perc. (1.7 10^-16) ; 40^th perc. (6.1 10^-9); 80^th perc. (4.0 10^-7); 95^th perc. (2.5 10^-6); 98^th perc. (5.0 10^-6); 100^th perc. (5.3 10^-5)

Fig. 8

Frame A illustrates the lava flow susceptibility at Nyiragongo volcano. The outlines of the 2002 (black line) and 2021 (hashed rectangles) lava flows of Nyiragongo volcano are superimposed on the lava flow susceptibility map in frame B. Frame C provides a zoom for the densely populated area of Goma and Gisenyi (black lines). Class boundaries are defined by susceptibility percentiles. Absolute susceptibility values for each class boundary are as follows: 0^th perc. (2.4 10^-16); 40^th perc. (7.1 10^-9); 80^th perc. (4.1 10^-7); 95^th perc. (3.8 10^-6); 98^th perc. (1.0 10^-5); 100^th perc. (2.5 10^-4)

At Nyiragongo, high to very high probabilities of lava flow invasion are found all around the base of the volcano, although to a lower extent on the northern and eastern flanks (Fig. 8A). The southern flank, including the areas invaded by lava during the 2002 and 2021 eruptions, is characterised by a high probability of lava flow invasion including the historical city centre of Goma (Fig. 8B). The map reveals that multiple channels with a very high probability of lava flow invasion extend towards the eastern part of Goma, with highest probabilities along major channels crossing the urban areas of Rukoko, Munigi and the zones destroyed by the main lava flow branches in 2002 and 2021. The probability of lava flow invasion is also very high in the western part of Goma with the presence of two main channels that can reach and invade this area. The central part of Goma, and Gisenyi in Rwanda, are less exposed to lava flow invasion (Fig. 8C).

The combination of the lava flow susceptibility of Nyamulagira and Nyiragongo volcanoes according to the two frequency ratio scenarios of Equation 2 is shown in Fig. 9. The first scenario assumes that the temporal probability of lava flow eruptions is the same at Nyamulagira and Nyiragongo volcanoes (W_Nya = W_Nyi = 0.5). It highlights that the probability of lava flow invasion increases by a few percentiles in the areas that can be invaded by lava flows originating from both volcanoes. This is clearly visible at the edge of the Nyamulagira and Nyiragongo volcanic fields, and within the inhabited areas located in the western part of Goma (Fig 9C and D). The lava flow susceptibility of areas that are only subjected to lava flow invasion by Nyamulagira or Nyiragongo volcano is only minorly affected by this integration

Fig. 9

Lava flow susceptibility in the western part of the VVP after combining the lava flow susceptibility maps of Nyamulagira and Nyiragongo volcanoes according to the two scenarios (scenario 1: W_Nya = W_Nyi = 0.5; scenario 2: W_Nya = 0.937 and W_Nyi = 0.063). Frame A and E show the resulting lava flow susceptibility map. Frame B and F show the difference in lava flow susceptibility percentile value compared to the maximum percentile value of the separate lava flow susceptibility maps of Nyamulagira and Nyiragongo volcanoes. Frame C, D, G and H provide a zoom on the area between the eastern part of Goma and Sake to highlight the changes in lava flow susceptibility in this inhabited area. Class boundaries are defined by susceptibility percentiles. Absolute susceptibility values for each class boundary are as follows: (A;C) 0^th perc. (8.3 10^-17) ; 40^th perc. (7.0 10^-9); 80^th perc. (2.8 10^-7); 95^th perc. (1.8 10^-6); 98^th perc. (3.9 10^-6); 100^th perc. (1.2 10^-4); (E; G) 0^th perc. (1.6 10^-17) ; 40^th perc. (5.3 10^-9); 80^th perc. (2.9 10^-7); 95^th perc. (2.0 10^-6); 98^th perc. (4.2 10^-6); 100^th perc. (5.0 10^-5)

The second scenario assumes that eruptions at Nyamulagira volcano are 15 times more likely than eruptions at Nyiragongo volcano. Although much more importance is provided to eruptions from Nyamulagira volcano when combining both susceptibility maps, the main channels going across the urban area of Goma remain at high to very high invasion probabilities. In this scenario the probability of lava flow invasion in the inhabited areas located in the western part of Goma decreases slightly compared to the individual lava flow susceptibility map (Fig. 9G and H).

Discussion

Lava flow susceptibility along Lake Kivu’s shoreline

Similarly to the study of Favalli et al. (2009), the lava flow susceptibility map of Nyiragongo volcano (Fig. 8) indicates that the areas invaded by the 2002 and 2021 lava flow eruptions - on the southern flank of Nyiragongo volcano - are characterised by a high to very high probability of lava flow inundation (Figs. 1 and 8C). The presence of the 2002 lava flows in the TanDEM-X DEM does not significantly modify the path of the lava flows, compared to simulations conducted on the 2000 SRTM DEM, due to the thin lava flow thickness associated with the low viscosity of Nyiragongo lavas (Morrison et al. 2020). The validity of this pattern was painfully confirmed by the emplacement of the two lava flow lobes of the May 22^nd 2021 eruption, that match with the zones with highest probabilities. These areas will likely continue to be characterised by high probabilities of lava flow invasion in the coming decades, even after the occurrence of multiple lava flows.

Our map also highlights high probability of lava flow invasion in the western part of Goma, where two main topographic channels could be invaded by lava to flow across the urban area (Fig. 8C). The latter was not observed in the lava flow susceptibility map of Favalli et al. (2009), as they restricted the occurrence of eruptive vents to the foot of the Nyiragongo main edifice and the north-south fracture system. In our study, the vent opening susceptibility at Nyiragongo was constrained by considering all the eruptive fissures and vents present in the volcanic field, without providing more weight to historical eruptive fissures and vents. The resulting vent opening susceptibility map indicates a secondary maximum of vent opening susceptibility in the south-western part of the Nyiragongo volcanic field, along the Rushayo alignment (Fig. 4B) and explains - together with the local topography - the increased lava flow inundation susceptibility in the western part of Goma. Lava flows of Nyamulagira volcano have also the potential to reach and invade the western part of Goma, albeit with a low to intermediate probability (Fig. 7A). As the classes of the susceptibility maps are defined in such a way to concentrate on the highest susceptibility values, one should consider that low and intermediate classes still represent areas that can potentially be impacted by lava flow in the future.

The urban area of Goma is densifying and sprawling not only to the north, towards Nyiragongo volcano, but also towards the west, along the northern shoreline of Lake Kivu (Michellier et al. 2013; Pech et al. 2018). The western part of Goma is characterised by high levels of population vulnerability (Michellier et al. 2020) and has been considered since 2002 as a potential development zone for the city of Goma. The presence of a high lava flow inundation susceptibility in that zone should be considered carefully in contingency planning and for the future development and expansion of Goma. High probability of lava flow inundation at several points along the Goma-Sake road also stresses the need to consider road blockage and evacuation route redundancy in emergency scenarios.

Lava flow susceptibility from neighbouring volcanoes

The vent opening susceptibility map of Nyamulagira and Nyiragongo volcanoes shows a small but significant spatial overlap close to the Nyamulagira and Nyiragongo main edifices and in the south-western part of the Nyiragongo volcanic field (Fig. 4C). It highlights that eruptive fissures and vents of both volcanoes can occur within that area. This overlap might be partly attributed to the diffusive effect of the kernel analysis, especially on the upper flank of the respective volcanoes, i.e. it is unlikely that magma from Nyiragongo would erupt on the upper flank of Namulagira, and vice-versa, but overlap between both volcanic systems along the lower flanks is realistic. This is indeed confirmed by the overlapping pattern of eruptive fissures and vents that could be associated to Nyamulagira and Nyiragongo volcanoes (Barette et al. 2016) based on their contrasted geochemical characteristics (Rogers et al. 1998; Chakrabarti et al. 2009). Since a geochemical characterisation is not available for each volcanic vent and fissure, they were attributed to Nyamulagira or Nyiragongo volcano according to the existing geological map (Thonnard et al. 1965). A detailed geochemical characterisation of all eruptive fissures and vents would allow better association of them to Nyamulagira and Nyiragongo volcanoes, and therefore produce a more accurate vent opening susceptibility map. Small changes in attribution of the eruptive fissures and vents to Nyamulagira and Nyiragongo volcanoes would however likely not modify the observed lava flow probability pattern significantly, due to the large number of closely spaced vents and the limited differences in lava flow input parameters for the two volcanoes. The orientation of the main volcanic rift zones of both volcanoes, and the distribution and orientation of distal fissure systems have been shown to be the result of the interactions of the volcanic and regional tectonic stress fields (Wadge and Burt 2011; Smets et al. 2015), such that magma intrusions sourced from either of the two volcanic systems could feed eruptive vents in the saddle zone between both edifices.

The calibration of the input parameters of the Q-LavHa model (Table 1, Fig. 6) for Nyiragongo and Nyamulagira lava flows results in differences in the Hc parameter and flow length distribution. The shorter length and lower Hc values obtained for Nyiragongo lava flows are consistent with the lower volume, lower viscosity, lower thickness and narrower flow lobes typically observed at Nyiragongo compared to Nyamulagira (Smets et al. 2015; Morrison et al. 2020). Figures 7 and 8 show that lava flow inundation probability around both volcanoes displays similar radial patterns, with high probability concentrating in multiple topographically low channels. The extent of historical lava flows confirms that high probabilities are modelled in zones recently occupied by lava flows. The map however also highlights zones with significant susceptibility in areas not recently affected by lava flows, i.e. in the western part of Goma, where the high lava flow invasion susceptibility overlaps with high susceptibility to explosive, phreatomagmatic eruptions and diffuse CO₂ degassing near the lake shore (Poppe et al. 2016b; Wauthier et al. 2018).

Figure 9 illustrates the added value of considering both volcanoes in the same model. A proper combination of the lava flow susceptibility of Nyamulagira and Nyiragongo is however challenging because the long-term temporal frequency of lava flow eruption is not well constrained. The scarcely available geochronological records and the presence of semi-permanent lava lakes at both volcanoes makes it hard to constrain the long-term temporal probability of lava flow eruptions at each volcano. Therefore, two different scenarios were developed to combine the lava flow susceptibility of Nyamulagira and Nyiragongo volcanoes. In this reference scenario assuming equal eruption frequency for both volcanoes, the area at the western edge of the city of Goma, along the border with the Virunga National Park, displays an increased probability of lava flow inundation, relative to the maps considering the volcanoes separately, as this area can be affected by lava flows coming from both volcanoes. Such insight is valuable for hazard communication and long-term urban planning. In the second scenario, the eruption history of the last 70 years of the two volcanoes was used to constrain a relative frequency of eruption, leading to a higher weight for the susceptibility map of the more active Nyamulagira volcano. The 2021 eruption of Nyiragongo, however, illustrates how sensitive such approach is for volcanoes with very few documented eruptions, and highlights the need to consider eruption frequency in such regional approach of volcanic hazard.

Limitation of the approach

The proposed integrated lava flow susceptibility map is affected by uncertainties and assumptions made in constraining the model. For the vent opening susceptibility map, fissures and vents were mapped based on their topographic expression and their identification on the geological map. Vents were attributed to a specific volcanic system based on the geological map but whether they were of historical or pre-historical age was not considered. Inaccuracy in the identification of specific vents should not drastically affect the general spatial pattern of the vent opening probability. Moreover, in areas where eruptions occur frequently along the same or overlapping fissure systems, the density of vents might be underestimated due to preferential overprinting of the vents. This is most likely to occur along the main volcanic fissure systems on the volcanoes upper flanks. Further geochemical and geochronological characterization of volcanic cones and fissures along the Rushayo alignment would also be needed to further constrain the volcanic system involved in their activity and the frequency of eruption in that area with high potential risk for the population (Barette et al. 2016; Poppe et al. 2016b).

The lava flow model considers the difference in lava flow length between the two volcanoes, but no variation in flow length within each volcanic field. Smets et al. (2015) identified that eruption volumes at Nyamulagira are increasing with distance from the volcano summit, but this is not directly translated in longer lava flows as it is moderated by lower background slopes. The number of flow units from Nyiragongo is too limited to constrain a spatial variation in length. Considering the distribution of vents and the topographic constraints, modification of the lava flow length parameter would only influence the regional susceptibility map north of Nyamulagira, but not along the densely inhabited shoreline of Lake Kivu. For crisis management, a better understanding of the factors controlling lava flow length is however needed to constrain potentially impacted zones based on exact eruptive vent location.

Finally, the presented lava flow susceptibility map is produced based on lava flow simulations run over a 30 m spatial resolution DEM derived from TanDEM-X data, representative for the topography in 2012 (Albino et al. 2015). The spatial resolution of the DEM was defined to reach the best balance between the spatial details of the model output and the computation power required to run the probabilistic model over a large area. The output probability map also realistically represents the width of zones typically affected by lava flow lobes. Running the model at high spatial resolution at local scale, i.e. within an urban area, would be possible as long as a highly accurate, recent DEM would be available and the interactions between lava flows and man-made infrastructure (buildings, roads) are properly understood (Tsang et al. 2020).

Emplacement of new lava flows modifies the topography and might cause such lava flow susceptibility map to become obsolete. Updating such map when new constraints on lava flow parameters or new or high accuracy topographic datasets become available is advisable. The model output nevertheless shows that zones recently covered by lava flows remain characterized by high susceptibility, suggesting that these high probability zones are not majorly modified by single lava flow emplacement at volcanoes with extremely low-viscosity lavas like Nyiragongo. Sensitivity tests (Fig. 13 in Appendix 2) comparing model outputs when using the SRTM DEM, representing the topography in 2000, instead of the 2012 DEM similarly showed limited change in the probability pattern, generally limited to local offsets of probability maxima to the lateral edge of lava flows in their distal area, where their thickness exceeded 2-3 meter. The produced susceptibility maps integrate constraints from the eruption history of the last century. The evolution of the topography over that time period has been limited to the formation of few scoria cones and the emplacement of individual lava flows with limited overlap. It is therefore reasonable to consider that the produced map is valid at a regional scale for several decades to come. We call for caution, however, in using such a map at a local scale, where micro-topography and the influence of anthropic modifications (urban fabric) should be taken into account.

Implications for volcanic risk management in Goma city

As demonstrated by the impact of the 2002 eruption and the chaotic evacuations during the extended volcano-seismic crisis of May 2021, long term reduction of the volcanic risk and efficient crisis management planning are essential to guarantee the safety of the growing population of the city of Goma (Smittarello et al. 2022). The lava flow susceptibility map presented here was developed within the GeoRisCA project (2012-2017), together with maps of population density and population vulnerability, and combined into a lava flow risk map (Michellier et al. 2020). The final versions of the map (Appendix 3) were presented and extensively discussed in a workshop in 2018 with a group of representatives from the Civil Protection, Goma Volcano Observatory (GVO) scientists, relevant urban planning and crisis management stakeholders and the authorities of Goma city and North Kivu province. One key conclusion was that the lava flow susceptibility map – assuming equal eruption probability between both volcanoes –, overlaid with key landmarks and the limits of the city’s administrative districts, was identified as the most easily understandable “volcanic eruption risk map” for Goma and its surroundings. Compared to a standard lava flow risk map combining information about lava flow susceptibility with population density and vulnerability, the lava flow susceptibility map alone was considered less technical and more appropriate for dissemination among the population to raise volcanic risk awareness (CAFOD-MRAC workshop in Goma, GVO, November 2018).

Special attention was also given to the choice of colours for the final lava flow susceptibility map (Thompson et al. 2015). Several tests of coloured maps were presented to GVO colleagues and other stakeholders. The two main arguments for choosing a “yellow to orange” scale are: (1) it is close to the colour of the lava, (2) it highlights nuances in the affected areas (i.e. feeling that not the whole area is dangerous); more nuanced than if the scale used is, for instance as presented during the workshop, from light to dark red. Regarding the final layout of the lava flow susceptibility map, the Civil Protection, GVO scientists and other involved stakeholders requested that additional information be included to improve the comprehension of the map: (1) the administrative boundaries and names of the districts of Goma, (2) Lake Kivu, the international border with Rwanda and the Virunga National park limits, and (3) specific landmarks. They also asked for a north-south orientation of the map. Finally, they decided that the short text displayed on the map to explain its content be presented in French only (and not in Swahili), to avoid any linguistic confusion.

After this workshop, and thanks to a close collaboration between the concerned stakeholders (mainly GVO, Civil Protection, Geographic Institute of Congo, NGO Caritas International – CAFOD), the map was made publicly available: it was printed in large format and displayed in 60 secondary schools of Goma in February 2019. It was also planned by the Civil Protection to take into consideration these maps and integrate them in both the contingency plan and the community evacuation plan for a volcanic eruption in Goma and its surroundings, but this process was not completed when the May 2021 eruption struck the city. Since then, the dissemination phase has resumed with the organisation of a workshop with the authorities of the city of Goma city and the province of North Kivu to further explain the lava flow susceptibility map and provide a large-format printed copy on canvas (to ensure its longevity) for consultation by the population in each district and commune office. However, even if the map could have been disseminated widely before the 2021 eruption, it is but one element to contribute to a more efficient and pro-active mitigation and management of future eruption risks, with the communication structure between local authorities playing a key role (Lowenstern et al. 2022).

Conclusion

We present and test a method to assess lava flow susceptibility from multiple closely-spaced volcanic systems with contrasted eruption characteristics. Using the case study of Nyiragongo and Nyamulagira volcanoes in the Virunga Volcanic Province, we highlight that specific attention needs to be dedicated to 1) areas were eruptive vents from both volcanic systems can open; 2) the contrast in the geometrical properties (i.e. length, thickness, width) of the lava flows; and 3) the contrasted eruption frequency of the two volcanic systems.

The produced lava flow susceptibility map for the Nyiragongo-Nyamulagira volcanic region highlights multiple channels with high lava flow inundation susceptibility affecting the city of Goma and the Goma-Sake road along the shoreline of Lake Kivu. Most of these places correspond well with the paths of the lava flows emplaced during the three historical eruptions of Nyiragongo (1977, 2002, 2021) and the multiple historical lava flows of Nyamulagira volcano. The western districts of Goma city and the zone of its potential westward extension, which have not been affected by eruptions in the last century, appear to be exposed to a significant lava flow susceptibility from Nyiragongo, as well as from Nyamulagira.

Although the integration of the lava flow inundation susceptibility from both volcanoes does not significantly modify the spatial pattern of the hazards, it attracts attention to areas exposed to hazards from both volcanoes. Considering the relative frequency of eruption from both volcanic systems also puts a new perspective on the volcanic hazards represented by Nyiragongo. The large population exposed to Nyiragongo more thancompensates, in terms of risk, for the lower eruption frequency observed in the last 70 years. Further constraints on the geochronology of the pre-historical eruptions from volcanic vents in the Nyamulagira-Nyiragongo volcanic field, and geochemical data to allocate them to a specific volcanic system, would help to further constrain in time the lava flow hazard and the associated risk.

As dramatically demonstrated by the human impact of the 2002 and 2021 eruptions of Nyiragongo, future eruptions of Nyamulagira and Nyiragongo volcanoes could cause major casualties and damages to the urban environment, and aggravate the fragile social, economic and political situation resulting from more than 20 years of armed conflicts (Reyntjens 2009). The further development and expansion of the urban area of Goma towards the north, in the direction of the Nyiragongo volcano, and the west, along the northern shoreline of Lake Kivu (Michellier et al. 2013; Pech et al. 2018) needs to consider the long-term spatial distribution of lava flow susceptibility from Nyiragongo and Nyamulagira volcanoes. The map presented here should inform this urban planning and serve as essential input for revising the volcanic risk management, including the revision of evacuation plans.

Appendices

Appendix 1

The computational time of Q-LavHA is proportional to the number of lava flow paths and inversely proportional to the distance between the vents (see Mossoux et al. 2016). We hence conducted a sensitivity test of the model to determine the number of lava flow paths and the distance between the vents that provides an optimal balance between the accuracy of the lava flow simulation and the computational time. The sensitivity test was performed on an area located on the eastern flank of Nyiragongo volcano (Fig. 10).

Fig. 10

Shaded relief of Nyiragongo volcano showing in a rectangle the source area considered for vents in the sensitivity test

The sensitivity test consists of analysing the relative difference (RD) between each overlapping pixel i of a non-computationally intensive simulation (S) and a reference simulation (Sref) as

$${RD}_{i}=\left(\frac{{S}_{i}-{S}_{ref,i}}{{S}_{ref,i}}\right)\times 100$$

Only pixels that were invaded by lava in both simulations were considered since it is not possible to attribute an error to a pixel that was invaded by lava in only one of the two simulations.

The non-computationally intensive and reference simulations of the sensitivity test are illustrated in Fig. 11. The results of the sensitivity test are summarized in Fig. 12. The results indicate that when lava flow paths are simulated from each pixel, the output of a lava flow simulation with one hundred or five hundred lava flow paths is more or less equivalent.

Fig. 11

Illustration of the lava flow simulations that were simulated with a distance between the vents of 30m (i.e. every pixel), 90m (i.e. every 3 pixels) and 150 m (i.e. every 5 pixels) with 50, 100, 200 and 500 iterations (500 iterations only for a distance of 30 meter between the vents). The reference simulation consists of a simulation with a distance between the vents of 30 m and 500 iterations

Fig. 12

Relative difference between each simulation and the reference simulation. The small inset provides a detail on the relative difference of the simulations realised with a distance between the vents of 30 meter

Appendix 2

Fig. 13

Comparison of the Q-LavHa simulation from the lowermost fissure of the 2002 eruption carried out on the SRTM DEM acquired in 2000 (A) and the TanDEM-X DEM acquired in 2012 (B), both at 30m spatial resolution. The topography of the 2002 lava flow (outlined in black), is present in the TanDEM X DEM but not in the SRTM DEM. Although differences can be noted between the simulation in the central and especially in the most distal part of the flow, the general areas susceptible to lava flow invasion remain very similar

Appendix 3

Fig. 14

Susceptibility map of lava flow inundation for the region of Goma city for the next eruption of Nyiragongo. The figure displays the final map (in french) as it was distributed in Goma after consultation of local stakeholders