# Subsurface lateral magma propagation from Nyiragongo volcano in the Western Rift Zone of the East African Rift

Virginie Pinel, Catherine A Mériaux

► To cite this version:

Virginie Pinel, Catherine A Mériaux. Subsurface lateral magma propagation from Nyiragongo volcano in the Western Rift Zone of the East African Rift. Journal of African Earth Sciences, 2025, 226, pp.105569. 10.1016/j.jafrearsci.2025.105569. ird-04995415

HAL Id: ird-04995415

https://ird.hal.science/ird-04995415v1

Submitted on 18 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Subsurface lateral magma propagation from Nyiragongo volcano in the Western Rift Zone of the East African Rift

Virginia Pinel a,*, Catherine A. Mériaux b

## ARTICLE INFO

**Keywords:**

Lateral magma propagation

Rift zone

Modelling

## ABSTRACT

Lateral magma propagation is a common feature of rift zones, with opening against the minimum compressive stress and vertical dykes flowing parallel to the rift direction. Depending on the competition between vertical and lateral magma migration, these dykes may either feed an eruption or not. In this context, the topography which includes the edifice load acts against the rise of the magma and favours lateral migration radially away from the edifice central area, thus feeding peripheral vents. Here, we focus on the Nyiragongo volcano, a volcanic edifice located in the western branch of the East African Rift and culminating at 3,470 meters above sea level. In practice, we study the combined effect of the extension induced by the rifting of the EAR and the topographic loading of the Nyiragongo volcano on the orientation of the dyke propagation plane and on the balance between the lateral and vertical propagation of the magma in this propagation plane. Using analytical and numerical models taking into account the effect of topography and the rift-induced west-east extensional stress field, we show that the path of a dyke originating from the volcanic edifice is first influenced by the load of the volcano, leading to a radial propagation. Beyond 5 km, however, the rift-induced extensional stress field dominates leading to a north-south propagation towards Lake Kivu. These results are consistent with the path of the magma deduced from the position of the vents and geophysical observations for the last two eruptive events of the Nyiragongo volcano (2002 and 2021). On the other hand, the lateral propagation over more than 20 km, where the magma remains trapped beneath the lake, is shown to be controlled by depth-dependent extension and reduced magma buoyancy and, to a lesser extent, the downslope towards Lake Kivu and the slight increase in rift extension towards the south.

## 1. Introduction

In rift zones, (sub)vertical magmatic intrusions, commonly called dykes, typically propagate over long distances in a direction parallel to the rift while remaining trapped at depth. This is observed for example in the Hawaiian Islands (Fiske and Jackson, 1972), in the Canary Islands (Delcamp et al., 2012) or in Iceland (Einarsson and Brandsdottir, 1980; Sigmundsson et al., 2015, 2022; Parks et al., 2023; Sigmundsson et al., 2024). This behaviour has been inferred from both observational and theoretical studies (Rubin and Pollard, 1987; Behn et al., 2006; Buck et al., 2006; Hamling et al., 2009; Grandin et al., 2012; Urbani et al., 2017, e.g.). In the East African Rift, similar lateral magmatic intrusions at shallow crustal depth occurred during the last two eruptions of Nyiragongo volcano in 2002 and 2021 (Wauthier et al., 2012; Smittarello et al., 2022). The direct consequence in terms of risk assessment is that magma can reach the surface at a large distance from the volcanic edifice, impacting remote areas. In the case

of the Nyiragongo volcano, the possibility of magma being emitted into Lake Kivu, located approximately 17 km south of the volcano, needs to be considered, particularly in view of the limnic eruptions that could potentially be induced (Schmid et al., 2005).

The physical understanding of such lateral magma propagation and propagation trajectories relies on two principles (Anderson, 1951, pp 21–25). First, dykes open as a result of tensile loading exerted at the tip of the intrusion. Second, they select the direction of “least resistance” by opening along planes across which the confining pressure is least. Characterizing the least confining pressure is thus the key to predicting the dyke propagation plane and its orientation. In this plane, the direction of propagation then depends on the location along the edge of the dyke where fracturing in purely mode I (opening mode) can happen, i.e. where the stress intensity factor $K_I$ exceeds or equals the ambient fracture toughness $K_c$ (Lawn and W., 1975, pp 58). Analytic solutions for the stress intensity factor are available in simple cases. For

<footnote>ª Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTERre, 38000 Grenoble, France</footnote>

<footnote>ᵇ ICTP-East African Institute for Fundamental Research, University of Rwanda, Kigali, Rwanda</footnote>

<footnote>* Correspondence to: ISTERre Université Savoie Mont-Blanc, Campus Scientifique, 73376 Le Bourget du Lac Cedex.
E-mail address: Virginia.Pinel@univ-smb.fr (V. Pinel).</footnote>
example, in the case of a uniform normal stress (i.e pressure) $p_0$ exerted on a planar dyke of length $a$ in an infinite region, the stress intensity factor is given by $K_I = p_0\sqrt{a}$. More generally, however, normal stresses exerted on the dyke surface are non-uniform. The stress intensity factor then depends on the geometry and the stress distribution along the dyke, which can be decomposed into average normal stress and stress gradients. In the case of a two-dimensional circular intrusion, the contribution to the stress intensity factor from the average stress is the same in all directions. As a result, the direction of propagation is controlled by the direction of the maximum stress gradient (Tada et al., 2000).

In practice, these lateral dykes indicate on the one hand that the principal horizontal stress is equal to or less than the vertical stress, and therefore that the plane in which the confining pressure is the lowest and where propagation occurs is vertical, and on the other hand that in this plane, the stress intensity factor is greater than or equal to the ambient fracture toughness $K_c$ at the lateral leading edge of the dyke rather than the vertical leading edge. The plane of least compression, its orientation and the direction of propagation in this plane are thus controlled by the relative amplitude of the principal vertical and horizontal stresses. These stresses are a function of topography, magmatic overpressure, buoyancy and the extensional tectonic stress field and its dependence on depth, dyke geometry (shape, length) and the load exerted on the dyke edges by these stresses. As far as topography is concerned, the compressive stress induced by loading at shallow depth counterbalances the buoyancy of the magma and tends to inhibit its vertical propagation (Pinel and Jaupart, 2000). This favours lateral propagation towards the flank vents, whose lateral distance is expected to increase for lower rates of magma influx or lower viscosity (Mériaux and Jaupart, 1995; Pinel and Jaupart, 2004). Although the influence of rifting-induced extension, topographic loading and buoyancy considered separately is now well understood, the combined effects of the corresponding stresses on the dyke propagation plane and on the balance between lateral and vertical magma propagation in the propagation plane have however not been examined together, which motivated this study.

The aim of this paper is not to present a complete three-dimensional dynamic model of the dyke trajectory but rather to focus on the combined first-order effects of edifice loading, local topography, buoyancy and the regional extensional stress field on the dyke trajectory using the East African Rift and the Nyiragongo volcano and its latest eruption in 2021 as a test case. The key question of this study is to better understand what controls lateral propagation in a rift zone. Although this question has been addressed before (see for instance (Rubin, 1990)), this study is distinguished by its analysis of static stresses in 3 dimensions. More specifically, we will successively examine how rifting-induced extension and topographic loading define the dyke propagation plane and its orientation and how buoyancy and extension then condition lateral magma propagation rather than vertical propagation in the propagation plane.

We summarize here the main observations taken from the abundant literature on the East African Rift and the Nyiragongo volcano. In the East African Rift, volcanism migrated from north to south as the rift propagated south and split into different branches, the Eastern branch and the Western branch. The associated magmatism also migrated from southern Ethiopia/northern Turkana to northern Tanzania between 45 Ma and 8 Ma, as shown by George et al. (1998). Locally, it has also been observed that volcanism often shifted from the interior of the valley to the exterior of the valley as the latter developed into a graben (Maccaferri et al., 2014). Overall, the East African Rift exhibits numerous volcanic plateaus, chains and provinces (> 20). Despite being qualified as amagmatic, the Western branch features four volcanic regions, the Toro-Ankole region in Western Uganda, the Rungwe volcanic field in southwestern Tanzania, South Kivu in the Democratic Republic of Congo (DRC) and the Virunga volcanic chain distributed over Uganda, the Democratic Republic of Congo (DRC) and

Rwanda (Ebinger, 1989). Spreading rate along the 5000 km long East African rift decreases from 6.5 mm/year in the north of the rift to 1.5 mm/year in the south of the rift (Stamps et al., 2008). By computing vertically averaged deviatoric stresses arising from horizontal gradients of gravitational potential energy derived from the CRUST 2.0 density model, Stamps et al. (2010) estimated the range of deviatoric extensional stresses for the East African Rift between 5 MPa and 15 MPa, highlighting an increase from south to north. This study is consistent with the previous estimate obtained by Coblentz and Sandiford (1994) with a lithospheric density model accounting for long-wavelength geoid anomalies. In the Virunga Volcanic Province, the rift trend is around N20 °E and the spreading rate has an intermediate value of 2.8 mm/yr. 250 km further south, the spreading rate reaches 3.2 mm/yr, while 250 km north it is only 2.1 mm/yr (Stamps et al., 2008). This implies that the north-south lateral gradient of this spreading rate is around $2.2 \times 10^{-3} \text{ mm yr}^{-1} \text{ km}^{-1}$ in this area. Around the Virunga Volcanic Province, principally extensional stresses have been estimated to reach a maximum value around 10 MPa (see Figure 4 of Stamps et al. (2010)).

As parts of the Virunga volcanic chain, Nyamuragira and Nyiragongo volcanoes located in the east of the Democratic Republic of Congo (DRC) are among the most active volcanoes in Africa. Nyiragongo volcano is a strato-volcano. It lies 14 km southeast of the Nyamuragira shield volcano and culminates at 3470 m overlooking the city of Goma in the DRC, the city of Gisenyi in Rwanda and Lake Kivu shared by the two countries. It emits rare alkaline and strongly silica-undersaturated lavas (Platz et al., 2004) characterized by an extremely low viscosity (Giordano et al., 2007) and expected to be highly degassed through the connection to the persistent lava lake present in the summit crater (Sawyer et al., 2008; Walwer et al., 2023). In 2002 and 2021, the eruption of Nyiragongo volcano was associated with simultaneous drainage of the summit lava lake (Komorowski et al., 2002/2003; Walwer et al., 2023). It not only caused lava flows on the southern flank of the volcano, but also the propagation of a subsurface dyke along apparently the same north-south eruptive fissure network (Komorowski et al., 2002/2003; Smittarello et al., 2022). The fissures were either newly formed or had been reactivated since the previous eruption in 1977, characterized by the same pattern of north-south eruptive fissures (Demant et al., 1994; Durieux, 2002/2003). Geodetic and seismic studies could infer the depth, direction and size of the dyke which was also shown to reach the bottom of Lake Kivu. The characteristics determined for these two intrusions are very similar; the propagation was initially radial to the south-east before moving approximately north-south in line with the eruptive fissures and fault lineaments in the Kivu basin (Smets et al., 2016); the intrusion plunged from the volcanic edifice to a maximum of 10 km depth under Lake Kivu. The main differences are the identification of two dykes during the 2002 eruption compared to only one in 2021, and openings of a maximum of 1 m in 2002 compared to 2.5 m in 2021 (Wauthier et al., 2012; Smittarello et al., 2022). In 2021, mapping of lava flows and vents using satellite imagery Poland (2022) suggests that the magma propagated in a N160 °E direction for around 5 km before aligning in the rift direction (N20 °E) as indicated by the location of the seismicity and the inversion of the geodetic data (Wauthier et al., 2012; Smittarello et al., 2022). In addition to these eruptive events, several episodes of lateral migration of magma at depth (about 10 km below sea level) towards the northeast were evidenced by seismic swarms and significant synchronous falls of the lava lake (Barrière et al., 2022) between 2016 and 2020.

In this paper, we first present the modelling approach and its limitations in Section 2. We evaluate the stresses due to the volcanic edifice loading and the topography and compare them with the extensional stresses induced by rifting (Section 3). This exercise allows us to show that the stresses induced by the volcanic edifice loading and the topography indicate a maximum compression oriented vertically, thus a propagation in a vertical plane opening perpendicular to the minimum compression. The comparison between the stresses further
<figure><img src="image_3.png" /><figcaption>Fig. 1. Geometry and stress field notations considered (a) in the 2D axisymmetric analytical case and (b) in the fully 3D numerical model. When discussing the orientation of the vertical plane of magma propagation (Section 4), the 2D stress field is considered in the horizontal planes (blue planes). When discussing the direction of propagation in the vertical propagation plane, determined by the force ratio between horizontal and vertical stress gradient (Section 5), the 2D stress field in the vertical planes is considered (red planes). Note that for the conical edifice (in a), the special case of a radial propagation plane is shown but other planes can be considered using the analytically calculated stress field.</figcaption></figure>

indicates that the influence of the edifice is predominant up to 5 km from the top of the edifice and that beyond this distance the influence of the extensional stresses induced by rifting dominates. We then analyse the orientation of the vertical propagation plane as a function of the edifice loading or the topography and the extension of rifting (Section 4). The analysis of the principal stress trajectories shows the change in orientation of the vertical propagation plane from a radial direction from the edifice to an orientation aligned with the north-south direction of the rift where the extensional constraints induced by rifting dominate. Section 5 finally analyzes the preferential propagation direction in the vertical propagation plane. This analysis highlights the importance of the increase in extension with depth, which allows to maintain the lateral propagation of magma at depth rather than a rise of magma towards the surface. The results are in agreement with the course of the last two eruptions of Nyiragongo. Sections 6 discusses the results of this study and their limitations and Section 7 is the final conclusion.

## 2. Modelling approach

Since we do not model the three-dimensional trajectory of the dyke, our analysis is broken down into several nested steps. We first show that the propagation plane is vertical. We then decompose the stress analysis successively in the horizontal plane, perpendicular to the vertical propagation plane, to determine its orientation, then in the vertical plane given its orientation, to evaluate the direction of propagation of the magma in this plane. In doing so, we use both an analytical and numerical approach. In the analytical approach, the surface geometry and load are reduced to a cone representative of the Nyiragongo volcanic edifice. In the numerical approach, the surface geometry and load correspond to a edifice/lake topography. Fig. 1 details for each of these two approaches the geometry and components of the stresses in and perpendicular to the propagation plane, to which we will refer in the following sections. For the numerical models, we use the commercial software COMSOL (https://www.comsol.com/products).

All the stresses we consider are perturbations of the lithostatic and isotropic stress reference state defined by Mc Garr (1988), as

$$ \sigma_{xx} = \sigma_{yy} = \sigma_{zz} = \rho_c g(H_{ref} - z), $$

where $\rho_c$ is the crustal density, $z$ is the elevation above sea level and $H_{ref}$ is the reference elevation, $z = H_{ref}$ being the surface of reference at which lithostatic stress is zero (see Fig. 1). The stress perturbations we consider are induced by rift-induced extensional stresses, stresses resulting from loadings such as the Nyiragongo volcanic edifice, Lake Kivu or general topography, and stresses arising from buoyancy forces.

We adopt the convention that compressive stresses are positive and conversely tensile stresses are negative.

Besides the two combined 2D analyses, our approach involves some simplifications. The crust is modelled by a linearly elastic homogeneous medium. In doing so, we ignore nonlinear deformations and pre-existing fissures and faults. As shown by field observations (e.g. Galland et al. (2019)), inelastic effects can occur at the crack tip, but they are limited to shallow propagation and affect only a limited area around the dyke tip, so that elastic behaviour of the crust away from the damaged zone should still control magma propagation. Furthermore pre-existing faults appear to have second-order importance in controlling magma pathing, as their orientation relative to the stress field is optimized for shear motions, whereas magma emplacement requires opening (Rivalta et al., 2019). In the case of Nyiragongo volcano, there are evidences suggesting that preexisting faults play a role. In particular, Smittarello et al. (2022) proposed that the Nyabihu fault mapped by Smets et al. (2016) stopped the propagation of the dyke during the 2021 eruption. However, this cessation of propagation only occurred after the dyke had propagated laterally for 20 km to the south. It is therefore reasonable to assume that there was no major faulting effect over at least the first 20 km of the dyke's trajectory. Similarly, fissures were observed at the surface during the 2021 eruption, but they likely did not interfere with the 20 km subsurface lateral propagation of the dyke (Smittarello et al., 2022). The potential effects of crustal anisotropy and faulting are described in more detail in the discussion (Section 6).

The dykes are also assumed to be aligned with the maximum compressive stress, and therefore open against the minimum compressive stress, as established by numerous field and theoretical studies (Anderson, 1951; Koenig and Pollard, 1998; Roman and Jaupart, 2014; Davis et al., 2021). However dyke propagation out of the plane of the principal compressive stress can occur (Mériaux and Lister, 2002; Watanabe et al., 2002; Menand et al., 2010; Maccaferri et al., 2019). In this respect, our modelling approach is a simplification. It essentially results from the fact that our approach is not dynamic but static. A dynamic propagation model would be necessary to appreciate the deviation from the trajectory inferred in this paper, as we discuss further in Section 6. However we will show that the results of our modelling approach agree well with the dyke path derived from geodetic data (Wauthier et al., 2012; Smittarello et al., 2022).

## 3. Influence of Mount Nyiragongo's edifice on the local stress field

In this section, we first establish a conical edifice load approximation for Nyiragongo volcano and an edifice-lake topographic profile. We subsequently isolate the influence of the volcanic edifice on the local
stress field using the conical edifice load approximation and an analyt-
ical solution and then considering the edifice-lake actual topography
and a numerical model. The stresses induced by the conical volcanic
edifice and the topography indicate the plane of least compression
suitable for propagation. These stresses are compared to the magnitude
of the stresses induced by the rift.

**3.1. Nyiragongo-Lake Kivu topographic profiles**

We combine the Digital Elevation Model (DEM) obtained from the
Shuttle Radar Topography Mission (SRTM) downloaded from the USGS
website (https://earthexplorer.usgs.gov/) with the bathymetry of Lake
Kivu (Lahmeyer, 1998; Ross et al., 2014) to generate an elevation
map, E(x, y), above sea level in the surrounding of Nyiragongo volcano
including Lake Kivu and the lake water thickness, W(x, y), with a
spatial resolution of 100 m. Note that the lake surface is located at an
elevation of 1461 m. The local topography and shoreline of Lake Kivu
are shown in Fig. 2a.

As shown by elevation profiles displayed in Fig. 2(b), Nyiragongo
volcano crater is located around 3180 m above sea level, it has a radius
of approximately 500 m and it is surrounded by crater walls hanging
250 m above the crater. It can be noted that the local slope reaches
30° in the summit part of Nyiragongo volcano and remains greater
than 5° on the flank of the volcano until reaching a distance of 6–
7 km from the summit of the volcano, except in the north-northwest
direction where a plateau is reached at a distance of 5 km due to the
topography related to the neighbouring Nyamulagira volcano. Beyond
5 km, the slope decreases and depends on the orientation of the profile
considered. Overall, the topography keeps decreasing on the profiles
oriented N180 or N225 °E, i.e. on profiles oriented in the direction of
the rift towards Lake Kivu, whereas on the profiles in the same direction
but oriented north or north-east, the elevation decreases much less after
a distance of 5 km. By averaging all the profiles, it can be assumed
that the volcanic edifice of Nyiragongo volcano has a typical height
H<sub>e</sub> of 1450 m and a radial extent R<sub>e</sub> of 5 km and its shape is well
approximated by a cone lying at an elevation of 2000 m above sea level
and thus characterized by a mean slope of 16°.

**3.2. Stress perturbation induced by the surface loading**

**3.2.1 Conical edifice load approximation and analytic solution**

We first use an analytical solution in cylindrical coordinates for the
2D axisymmetrical case to estimate the stress perturbation induced by
the load of a conical edifice of height H<sub>e</sub>, radius R<sub>e</sub> and density ρ<sub>e</sub> lying
over an elastic half-space with Poisson ratio ν of 0.25 (see Appendix A),
with a free surface located at the reference elevation H<sub>ref</sub>. Note that in
the following we take the density of the edifice equal to that of the crust
(ρ<sub>e</sub> = ρ<sub>c</sub>). Fig. 3 shows the stress components at various depths below
the surface. In the vicinity of the axis the minimum compressive stress
is horizontal (σ<sub>rr</sub> ≈ σ<sub>θθ</sub> < σ<sub>zz</sub>) such that a vertical dyke is expected
to form. The minimum horizontal stress varies with the distance from
the axis. Moving away from the axis but remaining close to it (grey
area in Fig. 3), the minimum horizontal stress is the hoop stress σ<sub>θθ</sub>,
such that radial dykes are expected. At a greater distance from the
axis, the radial stress σ<sub>rr</sub> becomes the minimum horizontal stress so that
circumferential dikes should in principle form. However the difference
between σ<sub>θθ</sub> and σ<sub>rr</sub> is small which does not guarantee the transition
from radial to circumferential. The transition distance at which the
radial stress becomes less than the circumferential stress increases with
depth.

The stresses due to the edifice are compared with the range of
extensional stresses estimated between 6 MPa and 10 MPa around the
Virunga province (Stamps et al., 2010). At a distance greater than
1.5 to 2 times the edifice radius, the influence of the edifice is no
longer present. In comparison, rift-induced extensional stresses become
dominant over the horizontal stresses due to the conical edifice at

<figure><img src="image_4.png" /></figure>

Fig. 2. Topography around Nyiragongo volcano. (a) Outline of lake Kivu (shown in white at an elevation of 1461 m above sea level) superimposed on the shaded SRTM Digital Elevation Model. Colour lines are for the height profiles, shown in (b), starting at Nyiragongo volcano summit and oriented N0 and N180 (in grey), N45 °E and N225 °E (in blue), N90 °E and N270 °E (in green) and N135 °E and N315 °E (in red). The black line represents the path followed by the magmatic intrusion during the 2021 eruption as inferred from eruptive vents location (Poland, 2022) and by the location of the seismicity and the inversion of the geodetic data (Smittarello et al., 2022). The black box indicates the limits of the 80 × 80 km area considered in our numerical model. Coordinates are expressed in the UTM system (UTM zone 35S). (b) Elevation as a function of the distance from the volcano summit along the height profiles shown in (a). The black dotted line represents a conical edifice with a radius of 5 km and a height of 1450 m, located above a reference elevation H<sub>ref</sub> of 2000 m.

shallow depths (i.e. z/R<sub>e</sub> = 0.2) at a distance as short as 0.6 times the
radius of the edifice. At depths greater than z/R<sub>e</sub> = 0.5, the rift-induced
extensional stresses dominate the horizontal stresses regardless of the
distance.

**3.2.2. Actual topography and numerical model**

Subsequently, we use the finite element numerical model, COM-
SOL Multiphysics (https://www.comsol.com/products), to quantify the
local stress field perturbation, relative to the lithostatic stress field, in-
duced by the local topography, including the effects of the Nyamulagira
volcano edifice and the rift shoulder, and the lake water load. The
numerical domain considered is a box of dimensions 80 km x 80 km
x 80 km, centred at the point of geographic coordinates -1.677383N,
29.171688E (750287, 9831805 in UTM coordinates, see the extension
of this domain in Fig. 2a). The whole domain is considered as an elastic
and homogeneous medium with a Poisson's ratio ν of 0.25. The top
surface of the box is located at an elevation of H<sub>ref</sub> = 1461 m above
sea level while the depth extends 80 km below sea level. Note that
we take a vertical dimension much larger than the crustal thickness
(estimated by Chakrabarti et al. (2009) to be around 35 km) not to
<figure><img src="image_5.png" /><figcaption>Fig. 3. Stress perturbation induced by a conical load of height $H_c = 1450$ m, radius $R_c = 5$ km and density $\rho_c = \rho_c = 2700 \text{ kgm}^{-3}$ applied at the surface of reference $H_{ref} = 2000$ m of an elastic half-space of Poisson coefficient $\nu = 0.25$. The curves in solid black, dashed red and dotted blue represent respectively the vertical ($\sigma_{zz}$), radial ($\sigma_{rr}$) and hoop ($\sigma_{\theta\theta}$) stress in cylindrical coordinates. The range of extensional stresses estimated for the East African Rift between 6 MPa and 10 MPa is represented with horizontal dashed lines for comparison. (a) Stress at elevation $z = 1000$ m, which corresponds to a normalized depth of $z/R_c = 0.2$ (depth in respect of the reference elevation and normalized by the edifice radius of 5 km). (b) Stress at elevation $z = -500$ m, which corresponds to a normalized depth of $z/R_c = 0.5$. (c) Stress at elevation $z = -3000$ m, which corresponds to a normalized depth of $z/R_c = 1$. The grey area is where the minimum compressive stress is the hoop stress ($\sigma_{\theta\theta}$), such that vertical radial dykes are expected. The blue vertical line represents the location of the lake shoreline for a dyke oriented north-south.</figcaption></figure>

consider the lithospheric mantle but only to ensure that the boundary condition applied at the bottom of the box does not influence the calculated crustal stress field. The domain is discretized into 350,000 tetrahedral elements whose size is smaller (maximum element size of 200 m) around the area of interest including the summit of Nyiragongo and the path followed by the magma during the 2021 eruption (see Fig. 2). We impose a zero displacement condition at the bottom and a zero horizontal displacement condition on the lateral faces of the box. At the upper surface, a normal stress corresponding to the load of the medium of density $\rho_c = 2700 \text{ kgm}^{-3}$ is applied by considering the difference between the local elevation $E(x, y)$ and the elevation of reference $z = H_{ref}$.

$$ \sigma_{zz}(x, y, H_{ref}) = \rho_c g (E(x, y) - H_{ref}). \quad (1) $$

Note that when the local elevation is below the elevation of reference, an unloading is applied to account for the loading deficit compared to the reference state of stress. In the area corresponding to Lake Kivu, an additional pressure corresponding to the load of the water is applied to the surface and the vertical stress applied then becomes at $z = H_{ref}$:

$$ \sigma_{zz}(x, y, H_{ref}) = \rho_c g (E(x, y) - H_{ref}) + \rho_w g W(x, y), \quad (2) $$

where $\rho_w$ is the water density taken equal to $1000 \text{ kgm}^{-3}$.

Fig. 4 shows the stress $\sigma_n$ (in dotted blue) acting perpendicular to a vertical dyke propagating from the central area of Nyiragongo volcano either towards the south (left column) or following the path observed during the 2021 eruption as reported in Fig. 2a (right column). This stress, which acts against the opening of the dyke when positive, decreases when going from Nyiragongo central area towards Lake Kivu along both paths considered. In addition, when the magma reaches Lake Kivu, this decrease is accentuated, particularly at shallow depths. The decrease is larger along the observed path (the right column in comparison the left one) and at shallower levels (the top row in comparison with the bottom one). In practice, this decrease is of the order of 15 MPa by the time the dyke has propagated ~ 10 km (see Fig. 4b), which is of the same order of magnitude as the tectonic extensional stress. This is expected to promote lateral spreading as the stress acting against the dyke opening decreases as magma propagates towards the lake. The only exception to this decrease of $\sigma_n$ along the magma path is observed when the dyke changes in orientation from $N160^\circ$ E to $N20^\circ$ E at depth ($z = -3550$ m) but it should be noted that here only the topography is considered while the tectonic stress is ignored. As the maximum extension is oriented perpendicular to rift direction, an extension greater than 2 MPa is enough to counterbalance this small increase (see Fig. B.12 in Appendix B). Also note that the stress $\sigma_n$ perpendicular to the dyke remains less compressive than the other stresses, which is consistent with the orientation of the plane of propagation.

At the shallow level ($z = 450$ m, Fig. 4a,b) the components of the stress field are similar to those estimated for a conical edifice using the analytical solution. Similarly, the edifice loading stress is larger than the extensional stress field when the distance from the summit remains smaller than 5 km. However, the attenuation with depth of the topography-induced stress field perturbation is much less pronounced than in the case of a conical load, mainly due to the larger spatial extent of the load considered. It follows that at a depth of approximately 1 km below sea level, the stress field induced by its load remains of the same order of magnitude as the extensional stress, below the volcanic edifice.

To summarize, both the analytical calculation for an axisymmetrical conical volcanic edifice and the numerical calculation considering the more complex real topography show that magma should propagate radially (i.e. vertical plane of propagation) and the local stress field of the edifice dominates the regional tectonic stress in the immediate vicinity of the edifice and at shallow depths relative to the elevation of reference. From the analytical solution, laterally propagating magma from the central area is expected to move radially from the edifice axis. Lateral magma propagation from the edifice is favoured by a decrease
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Stress perturbation induced by the local topography and lake water load along two vertical surfaces initiating at Nyiragongo summit: the north-south vertical plane (grey dashed profile in Fig. 2a) for the left column (a, c, e) and the vertical surface following the magma path oriented N160°E during 5 km (grey area) and then N20°E (black line in Fig. 2a) for the right column (b, d, f). The solid black and dotted blue curves respectively represent the vertical ($\sigma_{zz}$) and normal ($\sigma_n$) stress, whereas the dashed red curve represents the horizontal stress component ($\sigma_{perp}$) orthogonal to $\sigma_n$. The range of extensional stresses estimated for the Virunga region between 6 MPa and 10 MPa is also shown with horizontal dashed lines for comparison. (a), (b) Stress at elevation $z = 450$ m, which corresponds to a dimensionless depth of 0.2 (depth in respect of the reference elevation and normalized by the edifice radius of 5 km), (c) Stress at elevation $z = -1050$ m, which corresponds to a dimensionless depth of 0.5, (d) Stress at elevation $z = -3550$ m, which corresponds to a dimensionless depth of 1. The blue vertical line represents the location of the lake shoreline and the black vertical line represents the change in orientation of the vertical plane when considering the observed trajectory of the magma.</figcaption></figure>

of the horizontal stress acting perpendicular to the dyke. Since the ori-
tation of a laterally propagating vertical dyke is strongly dependent
on the direction of the principal horizontal stress components, in the
following section we calculate the orientation of the horizontal stress
field at various depths by combining the effects of the volcanic edifice
loading and the extension of the rift.

4. Orientation of a laterally propagating vertical dyke in an exten-
sional stress field

We calculate the two-dimensional stress field in an horizontal plane
located at a given depth either using the analytical formula for a
conical load or the numerical model based on the Digital Elevation
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Stress field perturbation induced by the combined effect of varying tectonic extension oriented N110 °E, perpendicular to the rift direction, and the load of a conical edifice (radius $R_e$ = 5000 m, height $H_e$ = 1450 m and density $\rho_e$ = 2700 kgm$^{-3}$) lying above the reference level ($H_{ref}$) set at an elevation $z$ = 2000 m. Black lines are for the orientation of the maximum compressive stress inside the horizontal plane along which a vertical dyke fully controlled by the local stress field should align. The colour scale corresponds to the horizontal deviatoric stress, with higher values favouring stronger control of the stress field over the orientation of the dyke. The black cross indicates the location of the summit of Nyiragongo and the red line shows the magma path of the 2021 eruption. The left, middle and right columns represent the cases of no extension, 1 MPa extension and 3 MPa extension, respectively. The upper, middle and lower rows are calculated at altitudes of 1000 m, -500 m and -3000 m respectively.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Stress field perturbation induced by the combined effect of varying tectonic extension oriented N110 °E, perpendicular to the rift direction, and the topographic load applied from the SRTM numerical model and the bathymetry of Lake Kivu above the reference level ($H_{ref}$) set at an elevation $z$ = 1461 m. Black lines are for the orientation of the maximum compressive stress inside the horizontal plane along which a vertical dyke fully controlled by the local stress field should align. The colour scale corresponds to the horizontal deviatoric stress, with higher values favouring stronger control of the stress field over the orientation of the dyke. The black cross indicates the location of the summit of Nyiragongo and the red line shows the magma path of the 2021 eruption. The left, middle and right columns represent the cases of no extension, 1 MPa extension and 3 MPa extension, respectively. The upper, middle and lower rows are calculated at altitudes of 450 m, -1500 m and -3550 m respectively.</figcaption></figure>

Model (DEM) and bathymetric data set. The stress field components, $\sigma_{xx}$, $\sigma_{yy}$, $\sigma_{xy}$, are expressed in Cartesian coordinates $(x, y)$ such that $x$ is in the east-west direction and positive towards the east, and $y$ is in the north-south direction positive towards the north as shown in Fig. 1. We add a tectonic stress field characterized by an amplitude $\sigma_{ext}$ and oriented N110 °E, i.e. with an angle $\theta$ = -20° with respect to the x-axis (east-west direction) as illustrated in Fig. 1. We then solve for the principal components of the total stress field.

Fig. 5 shows the direction of the maximum horizontal compressive stress as well as the amplitude of the horizontal deviatoric stress

$$ \sigma_2 - \sigma_3 = \sigma_{max} - \sigma_{min} $$

as a function of cone edifice load and varying rift extension. The left column, without extension, illustrates the load effect of the edifice alone. Consistent with the stress component along a radial profile at a given depth (see Fig. 3 and its description), vertical dykes propagating laterally at depth from the central zone of the volcano should be radial near the central zone before becoming circumferential at some distance, with a more distant transition for deeper dykes. However, the forcing of reorientation by stress is greater at shallower depths, as shown by the greater deviatoric stress in the zone where the maximum compressive stress changes from radial to concentric. Extensional stresses due to rift-ing tend to reorient the direction of maximum stress along the direction of the rift (N20 °E) both near the top of the edifice when the amount

of extension increases and at some distance from the edifice where the stress field induced by its load is negligible as already described in Roman and Jaupart (2014). The stress field induced 2500 m below the reference level, at 500 m below sea level, considering an extension of 1 MPa, matches particularly well the path followed by the magma during the 2021 eruption, explaining a reorientation occurring at 5 km from Nyiragongo summit. This extension is smaller than the extension inferred by Stamps et al. (2010). According to our calculations, an extension of between 0.5 and 2 MPa, at a depth of between 0.3$R_e$ and 1 $R_e$, which corresponds to an elevation of between 500 m and -3000 m, produces a stress field consistent with the path of the magma when we consider the effect of loading a conical edifice with a radius of 5 km and a height of 1450 m.

Fig. 6 shows the same information as Fig. 5 but considering the combined effect of the local stress field due to the topography and water load of the lake, estimated numerically with the COMSOL model, and the rift extension. In this case, the effect of Nyiragongo edifice load appears to be mixed by other topographic effects induced by mass distribution above the reference level (here set at $z$ = 1461 m) over a larger spatial extent such that the radial pattern of the maximum horizontal compressive in the vicinity of the volcanic edifice is much less obvious even in absence of extension (left panel). In particular, the orientation of the magma path in the first 5 kilometers from the summit is consistent to the local stress field only when considering no
extension and a deep level ($z = -3550$ m, bottom left panel in Fig. 6). As expected the extension always induces a stress field consistent with the magma path at some distance from the edifice. In addition the area where we observe a change in magma direction is always characterized by maximum values of the horizontal deviatoric stress. Here, our calculations favour a deeper propagation of the magma, at a depth of between $0.5R_e$ and $R_e$, which corresponds to an altitude of between -1000 and -3550 m, and a slightly greater extension, between 2 and 3 MPa, to produce a stress field fairly consistent with the trajectory of the magma.

In this section, we studied the stress field distribution in horizontal planes at fixed depths and a given extension in order to discuss potential effects on the orientation of laterally propagating vertical dykes. However, the depth of propagation may change as the magma propagates while at the same time the extension is expected to vary both vertically (Grandin et al., 2012) and horizontally. In the next section, we address the issue of competition between the horizontal and vertical stress gradient within the propagation plane, which controls the evolution of the propagation front and hence its depth with a depth- and position-dependent extension along the rift.

## 5. Competition between horizontal and vertical stress gradient in the plane of propagation

In this section, we estimate the preferential propagation direction in the propagation plane as a function of the horizontal and vertical stress gradients following the approach proposed by Davis et al. (2021) and recently applied by Mantiloni et al. (2023) in the specific context of calderas. We thus consider a finite batch of magma, the shape of a circular penny crack of radius $c$ in a vertical plane, and search for the maximum of stress intensity factor $K_I$ at its edges. The stress intensity factor $K_I$ at a given position of a penny-shaped crack edge and subjected to a linear pressure gradient is given by Tada et al. (2000):

$$K_I = \frac{4}{3\pi} \Delta\gamma c \sqrt{\pi c} \quad (3)$$

where $\Delta\gamma$ is the linear pressure gradient along the diameter passing through the position. $K_I$ is thus maximum at the position corresponding to the maximum pressure gradient.

The pressure gradient inside the crack results from the difference between the magmatic pressure and the external stress field acting perpendicular to the vertical plane. The magmatic pressure results from the static pressure gradient plus the viscous pressure induced by the flow of magma. In the case of a static crack, the magmatic pressure is only controlled by the magma buoyancy. The external stress field depends on the lithostatic stress field and any deviatoric stress fields such as the tectonic stress field, stresses induced by the propagating dyke, local topography and any load applied to the surface (e.g. lake water load). It follows that the pressure gradient $\partial P/\partial l$ along a given direction $l$ can be expressed as the buoyancy gradient $\partial P_{buoy}/\partial l$ and the external pressure gradient $\partial\sigma_n/\partial l$, where $\sigma_n$ is the disturbance of the normal stress acting at the wall of the magmatic intrusion with respect to the lithostatic stress. As a result,

$$\frac{\partial P}{\partial l} = \frac{\partial P_{buoy}}{\partial l} - \frac{\partial \sigma_n}{\partial l}. \quad (4)$$

We calculate the pressure gradient in the vertical propagation plane, using a 2D Cartesian coordinate system with $s$ axis corresponding to the horizontal axis of the plane oriented positively away from the volcanic edifice and $z$ the vertical axis oriented positively upward (see Fig. 1). The buoyancy gradient is maximum in the vertical direction and null in the horizontal with the components of the gradient be given:

$$\frac{\partial P_{buoy}}{\partial s} = 0, \quad (5)$$

$$\frac{\partial P_{buoy}}{\partial z} = (\rho_c - \rho_m)g, \quad (6)$$

where $\rho_m$ is the magma density. Providing the magma is less dense than the surrounding crust, this gradient is positive, ensuring a larger value of the stress intensity factor on the upper edge of the magma crack and thus an upward propagation. The calculation of the buoyancy gradient is estimated locally on a grid characterized by a given spatial step, which allows changes in crustal density to be taken into account with depth.

The components of the external stress field $\partial\sigma_n/\partial s$ and $\partial\sigma_n/\partial z$ are estimated numerically from the stress field normal to the vertical propagation plane taking into account stress perturbations of the lithostatic stress field induced by topographic loading, lake water loading and the extensional stress varying along the rift and with depth. As previously considered, extension is oriented N110 °E, i.e. with an angle $\theta = -20°$ with respect to the x-axis (east-west direction). The amplitude of extension is assumed to increase linearly with depth, according to Byerlee's law (Brace and Kohlstedt, 1980) and the extensional stress component is defined as:

$$\sigma_{ext}(z) = -\delta(H_{ref} - z), \quad (7)$$

where $\delta$ is the linear vertical gradient corresponding to the limit of rock strength when the minimum principal stress $\sigma_3$ is horizontal. As explained in the introduction, the spreading rate recorded around Lake Kivu is 2.8 mm yr$^{-1}$ with a north-south lateral gradient of about $2.2 \times 10^{-6}$ mm yr$^{-1}$ m$^{-1}$ (Stamps et al., 2008). Such a lateral gradient in the extensional stress field indicates that the extension amplitude increases slightly towards the south. By scaling the lateral extension gradient with the spreading rate gradient, we therefore define the amplitude of the horizontal stress gradient, at each given elevation $z$, as $|\sigma_{ext}(z)| \times 2.2 \times 10^{-6}/2.8$, with $\sigma_{ext}(z)$ given by Eq. (7). Using a value of the vertical gradient of extension ($\delta$) between 1400 $Pam^{-1}$ and 2500 $Pam^{-1}$, this corresponds to extensional stresses of respectively 8.4 and 15 MPa at 6 km depth. We note that the horizontal gradient of extension applied is of the order of $10 \ Pam^{-1}$ i.e. two orders of magnitude less than the vertical gradient. We then add the contributions of buoyancy and the external stress field due to the topography and estimate the direction corresponding to the maximum stress gradient. For the buoyancy term, magma density at Nyiragongo is expected to vary from 2530 to 2650 $kgm^{-3}$ (Burgi et al., 2014) while the crustal density is not well constrained, so we consider density contrasts between crust and magma ($\rho_c - \rho_m$) ranging from 200 to 0 $kgm^{-3}$.

The results for the direction of the maximum pressure gradient and, consequently, the expected direction followed by the magma within the propagation plane, are presented in Fig. 7. Note that Fig. 7 shows results for the stress fields due the conical edifice and the topography given by, respectively, the analytical solutions presented in Appendix A and the numerical solutions presented in Section 3.2.2. The behaviour is the same in both cases, the difference being mainly due to the difference in the reference level $H_{ref}$ considered: for the analytical solution, this level is higher, which makes it possible to calculate the stress field at a shallower level, which is however not entirely realistic since in the vicinity of Lake Kivu, the altitude is close to 1461 m.

Below the volcanic edifice, the magma ascent is inhibited by the compressive stress field induced by the edifice load such that the magma is expected to propagate laterally at the effective level of neutral buoyancy, which is the maximum overpressure level as defined by Pinel and Jaupart (2004). The depth of this level depends on the balance between the external stress field vertical gradient and the density contrast between the magma and the surrounding crust. It follows that whatever the depth at which the magma was initially stored, either at depth or even directly in the lava lake, once the dyke is initiated, it will reach this effective level of neutral buoyancy. This results is consistent with the lava lake drop observed simultaneously with each intrusive event (Barrière et al., 2022).

As the magma propagates away from the volcanic edifice, the compressive stress field decreases in amplitude, leading to a rise of the magma. The magma moves towards the surface at a distance of
<figure><img src="image_9.png" /></figure>

**Fig. 7.** Influence of the magma pressure gradient on the direction of propagation inside the vertical propagation plane. The propagation plane considered is the one represented in Fig. 2 (a) (N160 °E for 5 km then N20 °E). Black arrows are for the direction of the maximum pressure gradient inside a penny shaped crack of radius 200 m and the colour scale is for the amplitude of the normal stress acting normal to the vertical dyke ($σ_v$). The vertical blue line is the position of Lake Kivu shoreline. Red lines represent the paths followed by the magma starting below the edifice at $z = 0$, $z = -1$, $z = -2$ and $z = -3$ km. (a) Result obtained considering a conical load (analytical solution) applied above a surface located at elevation $H_{ref} = 2000$ m, an homogeneous density contrast $ρ_c - ρ_m = 150 \text{ kgm}^{-3}$ with a vertical extension gradient $δ = 1400 \text{ Pa/m}$. (b) Same as (a) except that the density contrast is reduced to an homogeneous value of $100 \text{ kgm}^{-3}$ in the shallow crust (at 1500 below the reference level $z = H_{ref}$). (c) Results obtained considering the topographic and water loading (numerical solution) applied above a surface located at elevation $H_{ref} = 1461$ m, an homogeneous density contrast $ρ_c - ρ_m = 200 \text{ kgm}^{-3}$ with a vertical extension gradient $δ = 2500 \text{ Pa/m}$. (d) Same as (c) except that the density contrast is reduced to an homogeneous value of $100 \text{ kgm}^{-3}$ in the shallow crust (at 1500 below the reference level $z = H_{ref}$).

approximately 5 km from the centre of the volcano, which corresponds to the location of the vent feeding the lava flow in 2021 on the flank of the volcano.

Then, due to the reorientation of the dyke, which induces a decrease in the stress acting perpendicular to the propagation plane, the magma sinks back to depth in agreement with the results from the inversion of geodetic data (Wauthier et al., 2012; Smittarello et al., 2022). The subsequent lateral spreading of magma over more than 20 km is strongly

related to the increase in extension with depth. In fact, comparing Figs. 7 and C.13 of Appendix C, with similar density contrasts but with and without depth-dependent extension, respectively, shows that this is a key ingredient for extensive lateral spreading. To a lesser extent, an additional decrease in buoyancy at shallow levels is required as shown in panels b and d of Fig. 7. The smaller density contrast at shallow levels may be explained by smaller crustal densities beneath Lake Kivu as inferred from geophysical imagery (Tuluka, 2010). Also, the magma is
expected to be degassed in this open system and possible additional gas loss could have occurred during the lava flow emplacement leading to an increase in magma density. The density difference between degassed and ascending magma is considered between 70 and 120 kgm⁻³ (Burgi et al., 2014).

Note that when the magma reaches the lake shoreline lateral prop-
agation is further promoted by the steep bathymetry of the lake. This
can be seen in the stress field calculated numerically (panels c and
d of Fig. 7) with arrows for the maximum pressure gradient steeply
descending close to the shoreline. However, this effect is difficult to
clearly demonstrate because our approximation of the stress gradient
is derived with a spatial step of 400 m which does not allow us to
precisely follow stress changes at shallow depth.

6. Discussion

Computing the stress field resulting from the combined effect of topographic loading and tectonic extension applied on an homogeneous elastic medium, we have shown in Section 3 that the influence of the Nyiragongo volcanic edifice on the stresses is limited to distances from the summit of the order of magnitude of the edifice radius at shallow depth. Beyond this distance, the extensional stresses due to rifting would dominate the stresses induced by the edifice. In particular, considering the path followed by the magma as aligned with the maximum compressive stress, the propagation of the dykes should initially be radial over a radius of 5 km, the minimum horizontal stress being the circumferential stress σ₀₀, to then change direction to align with that of the extensional stresses. This is in fact entirely consistent with Figure 4 of Smittarello et al. (2022) in which the propagating dyke changed direction 5 km from the summit approximately at about 1 km depth from the base line of the edifice during the 2021 eruption. We note that the same behaviour was reported by Wauthier et al. (2012) in their figure 8a during the 2002 eruption. The control of tectonic extensional stresses on dyke opening is typical of rift zones (e.g. Wright et al., 2012).

However, the crust in the Kivu basin is far from being a homo-
geneous medium as numerous fractures have been mapped. Most of
them are strongly aligned with the direction of the rift, but some
have different directions (Wadge and Burt, 2011; Smets et al., 2016).
In addition, magma remaining trapped at depth, either from failed
eruptions or from the cooled part of ancient dykes that fed eruptive
vents, also produces crustal heterogeneities. According (Smets et al.,
2016) (see Fig. 7A), the eruptive fissures are characterized by two
main orientations, one aligned with the rift direction (around N20 °E)
and the other with the N160 °E direction, thus corresponding exactly
to the two successive orientations followed by the magma during the
2021 eruption. All these crustal heterogeneities are expected to induce
crustal anisotropy in agreement with the results of seismic imagery in
rift zones (Keir et al., 2011).

Inherited faults and dykes can potentially modify the stress field. This effect could be studied by numerical modelling. For instance, the interaction between propagating magma and faults can be mod-
elled and has the potential to stop magma propagation (Maccaferri et al., 2016), supporting the potential role of the Nyabihu fault in stopping magma propagation during the 2021 eruption, as proposed by Smittarello et al. (2022). In general however, the lack of geometric constraints on the structures at depth makes the modelling difficult. This being said, the role of faults should only become significant when the control of the local stress field is weak (low deviatoric stress) or does not favour a particular direction of propagation. In particular, the N160 °E orientation of the magma trajectory over the first five kilo- metres for the 2021 event cannot be explained either by the analytical calculation of the edifice load, since in this axisymmetrical case any radial direction is identically favoured, nor very easily by the numerical calculation considering the topographic load. The N160 °E direction actually corresponds to the alignment of Nyamuragira and Nyiragongo

volcanoes and the expected orientation of transfer faults in the Virunga Volcanic Zone as described by Wadge and Burt (2011). As most of the faults and structures are oriented in the direction of the rift, they will only reinforce the role of extensional tectonic stress by favouring the N20 °E orientation of the magma once it has escaped the stress field of the Nyiragongo volcanic edifice. The presence of old fractures or fissures near the vent opened in 2021 may have favoured the arrival of magma at the surface as it was rising to the base of the volcano, as predicted by our model (see Fig. 7). However, the equilibrium between the stress field and buoyancy favoured the forward lateral propagation of the magma, so only a small amount of magma was emitted while most of it continued to propagate laterally at depth towards the south.

The ability of the magma to align its trajectory along the direction of the maximum compressive stress and its potential sensitivity to crustal heterogeneities depend on its driving pressure. In fact, many studies (Mériaux and Lister, 2002; Watanabe et al., 2002; Menand et al., 2010; Maccaferri et al., 2019) have shown that a dyke that is initially poorly aligned with respect to the principal stress trajectory will eventually reorient itself depending on the balance between the deviatoric stress and the magma driving pressure. If the driving pressure of the magma is high relative to the deviatoric stress, the magma will take longer to realign with the maximum compressive stress. In practice, the trajectory of the magma will be reoriented at a greater distance. For lateral propagation, the main pressure exerted on the magma results from its initial overpressure in the magma reservoir. In the particular case of the Nyiragongo volcano, an upper limit of this magmatic over-pressure can be estimated from the level of the lava lake. If we consider that the 2002 eruption reduced the magmatic overpressure inside the storage zone beneath Nyiragongo to zero overpressure (which is an upper limit because this reservoir may have been depressurized by the eruption), and that the level of the lava then rose by 300 m before the 2021 eruption (Barrière et al., 2022), a calculation of the hydrostatic pressure leads to an initial magmatic overpressure of less than 8 MPa. This upper limit of 8 MPa for the magma overpressure is comparable to the magnitude of the extension, so a delay in reorientation could occur. This means that the extension amplitude values derived in Section 4 could be underestimated, leading to a better agreement with the values previously provided by Stamps et al. (2010). Only a dynamic propagation model could resolve the exact trajectory, taking into account the pressure gradient of the magma resulting from its driving pressure and the viscous pressure loss induced by the flow. This model is beyond the scope of this study, as there is currently no comprehensive 3D model that takes into account magma flow and crustal fracture and can resolve both the orientation of the propagation plane and the direction of propagation within that plane. However, the validity of the approximation we used in Section 5 based on the direction of propagation of a penny-shaped crack could be tested using a model resolving the propagation of magma in a given plane (e.g. Zia and Lecampion (2020)) or by carrying out analog experiments.

From historical eruptions at Nyiragongo, it appears that propagation
over large distances (more that 20 km) only occurs towards the south
whereas some vents were observed in the north as for instance in
1977 (Durieux, 2002/2003; Wadge and Burt, 2011) and the intrusions
recorded by seismic swarms were directed towards the north (Barrière
et al., 2022). Lateral propagation in both directions from the feeding
centre is a common feature in rift zones. It was recently observed in
Iceland (Sigmundsson et al., 2024). At Nyiragongo volcano, even if
both directions can be taken, the propagation towards the south seems
more favoured. This could be explained both by the topography, which
decreases more markedly towards the south due in particular to the
proximity of the Nyamulagira edifice to the north-west, and by the
slight increase in extension towards the south in this particular zone
of the rift.

Kervyn et al. (2024) used historical vents and lava flows com-
bined with lava flow simulations to derive a lava flow probabilistic
susceptibility map of the Nyamuragira-Nyiragongo volcanic complex,
which is a very useful tool for risk mitigation. In their approach, the spatial distribution of the probability of a vent opening during the next eruption was constrained on the basis of the mapping of historical and prehistoric eruptive vents and fissures. Despite its obvious interest, this strategy is limited by the small number of well-known past eruptions and ignores the physical processes controlling the location of vents, making it impossible to take into account the evolution of the stress field potentially resulting from past activity. As already pointed out by Rivalta et al. (2019), a physics-based model such as the one proposed here has the advantage of allowing a better prediction of vent locations once the model parameters have been constrained on the basis of past eruptions. Such a model could certainly improve probabilistic maps of lava flow susceptibility.

## 7. Conclusion

We show that a static analysis of the stress field resulting from the combination of the topographic and lake load and the extensional stress field in the southern west branch of the East African rift can explain the two main characteristics of magma propagation from the Nyiragongo volcano. First, the propagation of magma from the volcanic edifice is influenced by the topographic load of the volcano, leading to a radial propagation over about 5 km beyond which the extensional stress field dominates leading to a north-south propagation towards Lake Kivu. Second, we show that magma subsurface propagation over more than 20 km requires a depth-dependent extension and a reduced density contrast at shallow depth, implying limited buoyancy, which is consistent with degassed magma propagating through layers of low-density upper crust together. Taking into account the influence of the stress field could potentially improve the assessment of the risks associated with the emplacement of lava flows. Finally, a similar subsurface propagation of magma towards Lake Kivu during future eruptions is likely and highlights that a subaquatic eruption at Lake Kivu with unknown consequences cannot be ruled out.

## CRediT authorship contribution statement

**Virginie Pinel:** Writing – original draft, Validation, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization. **Catherine A. Mériaux:** Writing – original draft, Validation, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

The present work was supported by the UNESCO International Geoscience Programme, IGCP project 767, EDOORs (https://e-doors.org, https://www.unesco.org/en/iggp/igcp-projects/767). We thank the two anonymous reviewers for their comments, which helped to improve the manuscript.

<figure><img src="image_10.png" /><figcaption>Fig. A.8. Stress Component $\sigma_{zz}$ due to a conical edifice in dimensionless form.</figcaption></figure>

## Appendix A. Stresses due to a conical edifice

The general form of the stresses induced by an axisymmetrical edifice of radius $R_e$ based on a surface $z = 0$ in a semi-infinite space can be found in Sneddon (1995) (p 469–473). Due to the axial symmetry, the stresses are expressed in cylindrical coordinates $(r, \theta, z)$, where the $z$ axis taken along the axis of symmetry, normal to the plane $z = 0$ and positive downwards. The boundary conditions are given in terms of the stresses $\sigma_{ij}$

$$ \sigma_{zz} = p(r) \text{ at } z = 0 \text{ and } r \le R_e, \quad (A.1) $$

$$ \sigma_{zz} = 0 \text{ at } z = 0 \text{ and } r > R_e, \quad (A.2) $$

$$ \sigma_{rz} = 0 \text{ at } z = 0, \quad (A.3) $$

$$ \sigma_{ij} = 0 \text{ as } z \to \infty, \quad (A.4) $$

where $p(r)$ is the pressure due to the load of the edifice. We note that we use the sign convention, positive compression and negative tension.

In an axially symmetric elastic medium unbounded in the $r$ direction, stresses and displacement can be expressed as a function of an arbitrary function $\Phi$, solution to the biharmonic equation $(\nabla^4 \Phi = 0)$. Stresses and displacements can then be given in the form of Hankel transforms (i.e. integrals of Bessel functions of the first kind $J_n(\xi r)$ of order $n$) of the function $\Phi$ satisfying the boundary conditions. In doing so, stresses $\sigma_{ij}$ and their trace $Tr(\sigma)$ can be expressed as follows:

$$ \sigma_{zz}(r, z) = z \int_{0}^{\infty} Z(\xi) J_{0}(\xi r) e^{-\xi z} \xi d\xi + \int_{0}^{\infty} Z(\xi) J_{1}(\xi r) e^{-\xi z} d\xi, \quad (A.5) $$

$$ \sigma_{rz}(r, z) = z \int_{0}^{\infty} Z(\xi) J_{1}(\xi r) e^{-\xi z} \xi d\xi, \quad (A.6) $$
<figure><img src="image_11.png" /><figcaption>Fig. A.9. Stress Component $\sigma_{rz}$ due to a conical edifice in dimensionless form.</figcaption></figure>

$$Tr(\sigma) = \sigma_{rr}(r, z) + \sigma_{\theta\theta}(r, z) + \sigma_{zz}(r, z) = 2(1 + \nu) \int_{0}^{\infty} Z(\xi) J_{0}(\xi r) e^{-\xi z} d\xi, \qquad (A.7)$$

$$\sigma_{rr}(r, z) = \int_{0}^{\infty} (1 - \xi z) Z(\xi) J_{0}(\xi r) e^{-\xi z} d\xi + \frac{1}{r} \int_{0}^{\infty} \left[ z - \frac{(1 - 2\nu)}{\xi} \right] Z(\xi) J_{1}(\xi r) e^{-\xi z} d\xi, \qquad (A.8)$$

where the function $Z(\xi)$ given by

$$Z(\xi) = \xi \int_{0}^{\infty} r p(r) J_{0}(\xi r) dr. \qquad (A.9)$$

has been introduced and $\nu$ is the Poisson's ratio. Full details of the derivation can be found in Sneddon (1995).

In the case of a conical edifice of density $\rho_e$ and height $h_e$, the pressure $p(r)$ is given by

$$p(r) = \rho_e g h_e \left(1 - \frac{r}{R_e}\right) \quad \text{for } r \le R_e, \qquad (A.10)$$

$$p(r) = 0 \quad \text{for } r > R_e, \qquad (A.11)$$

where $g$ is the gravitational acceleration. The first term corresponds to a pressure associated with a cylindrical load of height $h_e$ and radius $R_e$, the second term is a correction for the conical shape.

Substituting Eqs. (A.10) and (A.11) into Eq. (A.9) leads to

$$Z(\xi) = \xi \int_{0}^{R_e} r \rho_e g h_e \left(1 - \frac{r}{R_e}\right) J_0(\xi r), \qquad (A.12)$$

which can be expanded or simplified as

$$Z(\xi) = \rho_e g h_e \left[ R_e J_1(\xi R_e) - \xi \frac{R_e^2}{3} F_2\left(\frac{3}{2}; 1, \frac{5}{2}; -\frac{1}{4} \xi^2 R_e^2 \right) \right], \qquad (A.13)$$

<figure><img src="image_12.png" /><figcaption>Fig. A.10. Stress Component $\sigma_{rr}$ due to a conical edifice in dimensionless form.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. A.11. Stress Component $\sigma_{\theta\theta}$ due to a conical edifice in dimensionless form.</figcaption></figure>
Elevation z = -3550 m

<figure><img src="image_14.png" /><figcaption>**Fig. B.12.** Stress perturbation induced by the local topography and lake water load combined with a rift extension of 2 MPa along the vertical surface following the magma path as recorded by geophysical observations oriented N160 °E during 5 km (grey area) and then N20 °E (black line in Fig. 2a). The solid black and dotted blue curves represent respectively for the vertical (σ<sub>zz</sub>) and normal (σ<sub>n</sub>) stress, whereas the dashed red curve shows the horizontal stress component (σ<sub>perp</sub>) orthogonal to σ<sub>n</sub>. Stress is estimated at elevation z = -3550 m, which corresponds to a normalized depth of 1. The blue vertical line represents the location of the lake shoreline and the black vertical line represents the change in orientation of the vertical plane when considering the observed trajectory of the magma.</figcaption></figure>

where $\mathcal{F}_2(\frac{3}{2}; 1, \frac{5}{2}; -\frac{1}{4}\xi^2 R_e^2)$ is a generalized hypergeometric function.

Replacing Z(ξ) given by Eq. (A.13) into Eqs. (A.5)–(A.8) gives the stresses due to a conical edifice, which we can decomposed into the contribution for a cylindrical edifice of height h<sub>e</sub> and radius R<sub>e</sub> (Pinel and Jaupart, 2000; Roman and Jaupart, 2014) and the correction for the conical shape (terms involving the hypergeometric function). Eqs. (A.5)–(A.8) then become:

$$ \frac{\sigma_{zz}(r, z)}{\rho_e g h_e} = z R_e \int_0^\infty J_0(\xi r) J_1(\xi R_e) e^{-\xi z} \xi d\xi + R_e \int_0^\infty J_0(\xi r) J_1(\xi R_e) e^{-\xi z} d\xi \\ - z \frac{R_e^2}{3} \int_0^\infty \mathcal{F}_2\left(\frac{3}{2}; 1, \frac{5}{2}; -\frac{1}{4}\xi^2 R_e^2\right) e^{-\xi z} J_0(\xi r) \xi^2 d\xi \\ - \frac{R_e^2}{3} \int_0^\infty \mathcal{F}_2\left(\frac{3}{2}; 1, \frac{5}{2}; -\frac{1}{4}\xi^2 R_e^2\right) e^{-\xi z} J_0(\xi r) \xi d\xi \qquad (A.14) $$

$$ \frac{\sigma_{rz}(r, z)}{\rho_e g h_e} = z R_e \int_0^\infty J_1(\xi R_e) J_1(\xi r) e^{-\xi z} \xi d\xi \\ - z \frac{R_e^2}{3} \int_0^\infty \mathcal{F}_2\left(\frac{3}{2}; 1, \frac{5}{2}; -\frac{1}{4}\xi^2 R_e^2\right) J_1(\xi r) e^{-\xi z} \xi^2 d\xi \qquad (A.15) $$

$$ \frac{Tr(\sigma)}{\rho_e g h_e} = 2(1 + \nu) R_e \int_0^\infty J_0(\xi r) J_1(\xi R_e) e^{-\xi z} d\xi \\ - 2(1 + \nu) \frac{R_e^2}{3} \int_0^\infty \mathcal{F}_2\left(\frac{3}{2}; 1, \frac{5}{2}; -\frac{1}{4}\xi^2 R_e^2\right) J_0(\xi r) e^{-\xi z} \xi d\xi \qquad (A.16) $$

$$ \frac{\sigma_{rr}(r, z)}{\rho_e g h_e} = R_e \int_0^\infty (1 - \xi z) J_0(\xi r) J_1(\xi R_e) e^{-\xi z} d\xi \\ + \frac{R_e}{r} \int_0^\infty \left[ z - \frac{1 - 2\nu}{\xi} \right] J_1(\xi r) J_1(\xi R_e) e^{-\xi z} d\xi \\ - \frac{R_e^2}{3} \int_0^\infty (1 - \xi z) \mathcal{F}_2\left(\frac{3}{2}; 1, \frac{5}{2}; -\frac{1}{4}\xi^2 R_e^2\right) J_0(\xi r) e^{-\xi z} \xi d\xi \\ - \frac{R_e^2}{3r} \int_0^\infty \left[ z - \frac{1 - 2\nu}{\xi} \right] \mathcal{F}_2\left(\frac{3}{2}; 1, \frac{5}{2}; -\frac{1}{4}\xi^2 R_e^2\right) J_1(\xi r) e^{-\xi z} \xi d\xi \qquad (A.17) $$

<figure><img src="image_15.png" /><figcaption>**Fig. C.13.** Influence of the magma pressure gradient on the direction of propagation inside the vertical propagation plane. The propagation plane considered is the one represented in Fig. 2 a) (N160 °E for 5 km then N20 °E). Black arrows are for the direction of the maximum pressure gradient inside a penny shaped crack of radius 200 m and the colour scale is for the amplitude of the normal stress acting normal to the vertical dyke (σ<sub>n</sub>). The vertical blue line is the position of Lake Kivu shoreline. Red lines represent the paths followed by the magma starting below the edifice at z = 0, z = -1, z = -2 and z = -3 km. (a) Result obtained considering a conical load (analytical solution) applied above a surface located at elevation H<sub>ref</sub> = 2000 m, an homogeneous density contrast ρ<sub>e</sub> - ρ<sub>m</sub> = 150 kgm<sup>-3</sup>. No vertical extension gradient is considered. (b) Same as (a) except that the density contrast is reduced to an homogeneous value of 100 kgm<sup>-3</sup> in the shallow crust (at 1500 below the reference level z = H<sub>ref</sub>). (c) Results obtained considering the topographic and water loading (numerical solution) applied above a surface located at elevation H<sub>ref</sub> = 1461 m, an homogeneous density contrast ρ<sub>e</sub> - ρ<sub>m</sub> = 200 kgm<sup>-3</sup>. No vertical extension gradient is considered. (d) Same as (c) except that the density contrast is reduced to an homogeneous value of 100 kgm<sup>-3</sup> in the shallow crust (at 1500 below the reference level z = H<sub>ref</sub>).</figcaption></figure>

Eqs. (A.14)–(A.17) can be integrated numerically. In this study, we used Matlab and a in-house routine. The component σ<sub>θθ</sub> can be calculated from (A.16), (A.14) and (A.17).

Figs. A.8–A.10, show the dimensionless stresses components as a function of the dimensionless radial distance r for three dimensionless depths z = 0.2, 0.5 & 1. The stresses have been made dimensionless using the pressure scale ρ<sub>e</sub>gh<sub>e</sub>, where ρ<sub>e</sub> and g are respectively the density of the edifice and the gravitational acceleration. The distance r and depths z have been made dimensionless using the length scale R<sub>e</sub>. The Poisson ratio ν was taken equal to 0.25. We verified that solutions at the axis (r = 0) for the cylindrical and conical edifice corresponds to the analytical solution provided by Pinel and Jaupart (2000) by linearization.

Corrections to all stress components from the cylindrical shape to the conical shape decrease by about half the magnitude of the stress
components due to a cylindrical edifice. At a dimensionless distance $r/R_e$ between 1.5 and 2, the influence of the edifice is no longer present.

## Appendix B. Stress perturbation induced by the local topography and lake water load combined with a rift extension of 2 MPa along the vertical surface following the magma path

In this appendix, we show the combined effect of local topography and tectonic extension on the stress field along the magma path. In particular, as the dyke reorients in the direction of the rift, the normal stress acting on its wall decreases due to the greater influence of tectonic extension.

## Appendix C. Magma pressure gradient and direction of propagation inside the vertical propagation plane without depth-dependent extension

In this appendix, we show the importance of the depth dependence of rift-induced extension. Fig. C.13 shows the magma path with a density similar to the model of Fig. 7 but without depth-dependent extension. Comparison of these two Figs. 7 and C.13 demonstrates that depth dependence is a key ingredient for extensive lateral propagation. Without it, lateral magma propagation is significantly reduced to just over 5 km versus over 20 km propagation when taking into account the increase in extension with depth, regardless of the loading considered (conical edifice or topography).

## Data availability

Data will be made available on request.

## References

Anderson, E.M., 1951. The Dynamics of Faulting and Dyke Formation with Applications To Britain, second ed. Oliver and Boyd, Edinburgh.

Barrière, J., d'Oreye, N., Smets, B., Oth, A., Delhaye, L., Subira, J., Mashagiro, N., Der- auw, D., Smittarello, D., Syavulisembo, A.M., Kervyn, F., 2022. Intra-crater eruption dynamics at nyiragongo (D.R. Congo), 2002–2021. J. Geophys. Res.: Solid Earth 127, http://dx.doi.org/10.1029/2021JB023858, e2021JB023858 2021JB023858.

Behn, M.D., Buck, W.R., Sacks, I.S., 2006. Topographic controls on dike injection in volcanic rift zones. Earth Planet. Sci. Lett. 246, 188–196.

Brace, W.F., Kohlstedt, D.L., 1980. Limits on lithospheric stress imposed by laboratory experiments. J. Geophys. Res.: Solid Earth 85, 6248–6252. http://dx.doi.org/10.1029/BJ085B11p06248.

Buck, W.R., Einarsson, P., Brandsdóttir, B., 2006. Tectonic stress and magma chamber size as controls on dike propagation: Constraints from the 1975–1984 Krafla rifting episode. J. Geophys. Res.: Solid Earth 111, B12404. http://dx.doi.org/10.1029/2005JB003879.

Burgi, P.Y., Darrah, T.H., Tedesco, D., Eymold, W.K., 2014. Dynamics of the mount nyiragongo lava lake. J. Geophys. Res.: Solid Earth 119, 4106–4122.

Chakrabarti, R., Basu, A.R., Santo, A.P., Tedesco, D., Vaselli, O., 2009. Isotopic and geochemical evidence for a heterogeneous mantle plume origin of the Virunga volcanics, Western rift, East African Rift system. Chem. Geol. 259, 273–289. http://dx.doi.org/10.1016/j.chemgeo.2008.11.010.

Coblentz, D.D., Sandiford, M., 1994. Tectonic stresses in the African plate: Constraints on the ambient lithospheric stress state. Geology 22, 831–834. http://dx.doi.org/10.1130/0091-7613(1994)022<0831:TSITAP>2.3.C02.

Davis, T., Bagnardi, M., Lundgren, P., Rivalta, E., 2021. Extreme curvature of shallow magma pathways controlled by competing stresses: insights from the 2018 Sierra Negra eruption. Geophys. Res. Lett. 48, e2021GL093038.

Delcamp, A., Troll, V.R., Van Wyk de Vries, B., Carracedo, J.C., Petronis, M.S., Pérez- Torrado, F.J., Deegan, F.M., 2012. Dykes and structures of the NE rift of Tenerife, Canary Islands: a record of stabilisation and destabilisation of ocean island rift zones. Bull. Volcanol. 74, 963–980.

Demant, A., Lestrade, P., Lubala, R.T., Kampunzu, A.B., Durieux, J., 1994. Volcanological and petrological evolution of Nyiragongo volcano, Virunga volcanic field, Zaire. Bull. Volcanol. 56, 47–61.

Durieux, J., 2002/2003. Nyiragongo: The january 10th 1977 eruption. Acta Volcanol. 14–15, 145–148.

Ebinger, C.J., 1989. Tectonic development of the western branch of the East African rift system. Geol. Soc. Am. Bull. 101, 885–903.

Einarsson, P., Brandsdóttir, B., 1980. Seismological evidence for lateral magma intrusion during the 1978 deflation of the Krafla Volcano in NE Iceland. J. Geophys. 47, 160–165.

Fiske, R., Jackson, E., 1972. Orientation and growth of Hawaiian volcanic rifts: the effect of regional structure and gravitational stresses. Proc. R. Soc. A 329, 299–326.

Galland, O., Spacapan, J.B., Rabbel, O., Mair, K., Soto, F.G., Eiken, T., Schiuma, M., Leanza, H.A., 2019. Structure, emplacement mechanism and magma-flow significance of igneous fingers – implications for sill emplacement in sedimentary basins. J. Struct. Geol. 124, 120–135. http://dx.doi.org/10.1016/j.jsg.2019.04.013.

George, R., Rogers, N., Kelley, S., 1998. Earliest magmatism in Ethiopia: evidence for two mantle plumes in one flood basalt province. Geology 26, 923–926.

Giordano, D., Polacci, M., Longo, A., Papale, P., Dingwell, D.B., Boschi, E., Kasereka, M., 2007. Thermo-rheological magma control on the impact of highly fluid lava flows at Mt. Nyiragongo. Geophys. Res. Lett. 34, http://dx.doi.org/10.1029/2006GL028459.

Grandin, R., Socquet, A., Doubre, C., Jacques, E., King, G.C.P., 2012. Elastic thickness control of lateral dyke intrusion at mid-ocean ridges. Earth Planet. Sci. Lett. 83–95. http://dx.doi.org/10.1016/j.epsl.2011.12.011.

Hamling, I.J., Ayele, A., Bennati, L., Calais, E., Ebinger, C.J., Keir, D., Lewi, E., Wright, T.J., Yirgu, G., 2009. Geodetic observations of the ongoing Dabbahu rifting episode: new dyke intrusions in 2006 and 2007. Geophys. J. Int. 178, 989–1003.

Keir, D., Belachew, M., Ebinger, C.J., Kendall, J.M., Hammond, J.O.S., Stuart, G.W., Ayele, A., Rowland, J.V., 2011. Mapping the evolving strain field during continental breakup from crustal anisotropy in the Afar Depression. Nat. Commun. 2, 285.

Kervyn, M., Barette, F., Poppe, S., Smets, B., Syavulisembo Muhindo, A., Kam- bale Makundi, J., Gnungzi Kahashi, Y., Kambale Ndagana, J., Mossoux, S., Kervyn, F., Michellier, C., 2024. Assessing lava flow susceptibility at neighbouring volcanoes: Nyamulagira and Nyiragongo volcanoes, Virunga Volcanic Province. J. Appl. Volcanol. 13, 5.

Koenig, E., Pollard, D.D., 1998. Mapping and modeling of radial fracture patterns on Venus. J. Geophys. Res.: Solid Earth 103, 15183–15202.

Komorowski, J., Tedesco, D., Kasereka, M., Allard, P., Papale, P., Vaselli, O., Durieux, J., Baxter, P., Halbwachs, M., Akumbe, M., Baluku, B., Briole, P., Ciraba, M., Dupin, J.C., Etoy, O., Garcin, D., Hamaguchi, H., Houlie, N., Kavotha, K.S., Lemarchand, A., Lockwood, J., N., L., Mavonga, G., de Michele, M., S., M., Mukambilwa, K., Munyololo, F., Newhall, C., J., R., Yalire, M., Wafula, M., 2002/2003. The 2002 eruption-The 2002 flank eruption of Nyiragongo Volcano (Democratic Republic of Congo): chronology, evidence for a tectonic rift trigger, and impact of lava flows on the city of Goma. Acta Volcanol. 14–15, 27–62.

Lahmeyer, Osae, 1998. Bathymetric survey of lake kivu. Final report. Technical Report, Republic of Rwanda, Ministry of Public Work, Directory of Energy and Hydrocarbons, Kigali.

Lawn, B.R., W., T.R., 1975. Fracture of Brittle Solids. Cambridge University Press.

Maccaferri, F., Rivalta, E., Keir, D., Acocella, V., 2014. Off-rift volcanism in rift zones determined by crustal unloading. Nat. Geosci. 7, 297–300.

Maccaferri, F., Rivalta, E., Passarelli, L., Aoki, Y., 2016. On the mechanisms governing dike arrest: Insight from the 2000 Miyakejima dike injection. Earth Planet. Sci. Lett. 434, 64–74. http://dx.doi.org/10.1016/j.epsl.2015.11.024.

Maccaferri, F., Smittarello, D., Pinel, V., Cayol, V., 2019. On the propagation path of magma-filled dikes and hydrofractures: The competition between external stress, internal pressure, and crack length. Geochem. Geophys. Geosyst. 20, 2064–2081. http://dx.doi.org/10.1029/2018GC007915.

Mantiloni, L., Rivalta, E., Davis, T., 2023. Mechanical modeling of pre-eruptive magma propagation scenarios at calderas. J. Geophys. Res.: Solid Earth 128, http://dx.doi.org/10.1029/2022JB025956, e2022JB025956.

Mc Garr, A., 1988. On the state of lithospheric stress in the absence of applied tectonic forces. J. Geophys. Res.: Solid Earth 93, 13, 609–13, 617.

Menand, T., Daniels, K.A., Benghiat, P., 2010. Dyke propagation and sill formation in a compressive tectonic environnement. J. Geophys. Res.: Solid Earth 115, B08201.

Mériaux, C., Jaupart, C., 1995. Simple fluid dynamic models of volcanic rift zones. Earth Planet. Sci. Lett. 136, 223–240.

Mériaux, C., Lister, J.R., 2002. Calculation of dike trajectories from volcanic centers. J. Geophys. Res.: Solid Earth 107, ETG-10.

Parks, M., Sigmundsson, F., Drouin, V., Hjartardóttir, Á.R., Geirsson, H., Hooper, A., Vogfjörd, K.S., Ófeigsson, B.G., Hreinsdóttir, S., Jensen, E.H., et al., 2023. Deformation, seismicity, and monitoring response preceding and during the 2022 Fagradalsfjall eruption, Iceland. Bull. Volcanol. 85, 60.

Pinel, V., Jaupart, C., 2000. The effect of edifice load on magma ascent beneath a volcano. Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 358, 1515–1532.

Pinel, V., Jaupart, C., 2004. Magma storage and horizontal dyke injection beneath a volcanic edifice. Earth Planet. Sci. Lett. 221, 245–262. http://dx.doi.org/10.1016/j.epsl.2004.01.031.

Platz, T., Foley, S.F., André, L., 2004. Low-pressure fractionation of the Nyiragongo volcanic rocks, Virunga Province, D.R. Congo. J. Volcanol. Geotherm. Res. 136, 269–295. http://dx.doi.org/10.1016/j.jvolgeores.2004.05.020.
Poland, M.P., 2022. Synthetic aperture radar volcanic flow maps (SAR VFMs): a simple method for rapid identification and mapping of volcanic mass flows. Bull. Volcanol. 84, 32.

Rivalta, E., Corbi, F., Passarelli, L., Acocella, V., Davis, T., Di Vito, M.A., 2019. Stress inversions to forecast magma pathways and eruptive vent location. Sci. Adv. 5, http://dx.doi.org/10.1126/sciadv.aau9784.

Roman, A., Jaupart, C., 2014. The impact of a volcanic edifice on intrusive and eruptive activity. Earth Planet. Sci. Lett. 408, 1–8.

Ross, K.A., Smets, B., De Batist, M., Hilbe, M., Schmid, M., Anselmetti, F.S., 2014. Lake-level rise in the late Pleistocene and active subaquatic volcanism since the Holocene in Lake Kivu, East African Rift. Geomorphology 221, 274–285. http://dx.doi.org/10.1016/j geomorph.2014.05.010.

Rubin, A.M., 1990. A comparison of rift-zone tectonics in Iceland and Hawaii. Bull. Volcanol. 52, 302–319.

Rubin, A., Pollard, D., 1987. Origins of blade-like dikes in volcanic rift zones. US Geological Survey Professional Paper, 1350, pp. 1449–1470.

Sawyer, G.M., Carn, S.A., Tsanev, V.I., Oppenheimer, C., Burton, M., 2008. Investigation into magma degassing at Nyiragongo volcano, Democratic Republic of the Congo. Geochem. Geophys. Geosyst. 9, http://dx.doi.org/10.1029/2007GC001829.

Schmid, M., Halbwachs, M., Wehrli, B., Wüest, A., 2005. Weak mixing in Lake Kivu: New insights indicate increasing risk of uncontrolled gas eruption. Geochem. Geophys. Geosyst. 6, http://dx.doi.org/10.1029/2004GC000892.

Sigmundsson, F., Hooper, A., Hreinsdottir, S., Vogfjörd, K.S., Ofeigsson, B.G., Heimisson, E.R., Dumont, S., Parks, M., Spaans, K., Gudmundsson, G.B., Drouin, V., Arnadottir, T., Jonsdottir, K., Gudmundsson, M.T., Hognadottir, T., Fridriksdottir, H.M., Hensch, M., Einarsson, P., Magnusson, E., Samsonov, S., Brandsdottir, B., White, R.S., Agustsdottir, T., Greenfield, T., Green, R.G., Hjartardottir, A.R., Pedersen, R., Bennett, R.A., Geirsson, H., La Femina, P.C., Bjornsson, H., Palsson, F., Sturkell, E., Bean, C.J., Mollhoff, M., Braiden, A.K., Eibl, E.P.S., 2015. Segmented lateral dyke growth in a rifting event at Bárðarbunga volcanic system, Iceland. Nature 517, 191–195.

Sigmundsson, F., Parks, M., Geirsson, H., Hooper, A., Drouin, V., Vogfjörd, K.S., Ófeigsson, B.G., Greiner, S.H., Yang, Y., Lanzi, C., et al., 2024. Fracturing and tectonic stress drive ultrarapid magma flow into dikes. Science 383, 1228–1235.

Sigmundsson, F., Parks, M., Hooper, A., Geirsson, H., Vogfjörd, K.S., Drouin, V., Ófeigsson, B.G., Hreinsdóttir, S., Hjaltadóttir, S., Jónsdóttir, K., Einarsson, P., Barsotti, S., Horálek, J., Ágútsdóttir, T., 2022. Deformation and seismicity decline before the 2021 fagradalsfjall eruption. Nature 609, 523–528.

Smets, B., Delvaux, D., Ross, K.A., Poppe, S., Kervyn, M., d'Oreye, N., Kervyn, F., 2016. The role of inherited crustal structures and magmatism in the development of rift segments: Insights from the Kivu basin, western branch of the East African Rift. Tectonophysics 683, 62–76.

Smittarello, D., Smets, B., Barrière, J., Michellier, C., Oth, A., Shreve, T., Grandin, R., Theys, N., Brenot, H., Cayol, V., Allard, P., Caudron, C., Chevrel, O., Darchambeau, F., de Buyl, P., Delhaye, L., Derauw, D., Ganci, G., Geirsson, H., Kamate Kaleghetso, E., Kambale Makundi, J., Kambale Nguomoja, I., Kasereka Mahinda, C., Kervyn, M., Kimanuka Ruriho, C., Le Mével, H., Molendijk, S., Namur, O., Poppe, S., Schmid, M., Subira, J., Wauthier, C., Yalire, M., d'Oreye, N., Kervyn, F., Syavulisembo Muhindo, A., 2022. Precursor-free eruption triggered by edifice rupture at Nyiragongo volcano. Nature 609, 83–88.

Sneddon, I.N., 1995. Fourier Transforms. Courier Corporation.

Stamps, D.S., Calais, E., Saria, E., Hartnady, C., Nocquet, J.M., Ebinger, C.J., Fernandes, R.M., 2008. A kinematic model for the East African Rift. Geophys. Res. Lett. 35, http://dx.doi.org/10.1029/2007GL032781.

Stamps, D., Flesch, L., Calais, E., 2010. Lithospheric buoyancy forces in Africa from a thin sheet approach. Int. J. Earth Sci. 99, 1525–1533.

Tada, H., Paris, P., Irwin, G., 2000. The Stress Analysis of Cracks Handbook (Ed.). The American Society of Mechanical Engineers.

Tuluka, G.M., 2010. Crustal structure beneath two seismic broadband stations revealed from teleseismic P-wave receiver function analysis in the Virunga volcanic area, Western Rift Valley of Africa. J. Afr. Earth Sci. 58, 820–828. http://dx.doi.org/10.1016/j.jafrearsci.2009.11.003, active Volcanism and Continental Rifting in Africa.

Urbani, S., Acocella, V., Rivalta, E., Corbi, F., 2017. Propagation and arrest of dikes under topography: Models applied to the 2014 Bardarbunga (Iceland) rifting event. Geophys. Res. Lett. 44, 6692–6701.

Wadge, G., Burt, L., 2011. Stress field control on eruption dynamics at a rift volcano: Nyamuragira, D. R. Congo. J. Volcanol. Geotherm. Res. 207, http://dx.doi.org/10.1016/j.volgeores.2011.06.012.

Walwer, D., Wauthier, C., Barrière, J., Smittarello, D., Smets, B., d'Oreye, N., 2023. Modeling the intermittent lava lake drops occurring between 2015 and 2021 at Nyiragongo volcano. Geophys. Res. Lett. 50, e2022GL102365.

Watanabe, T., Masuyama, T., Nagaoka, K., Tahara, T., 2002. Analog experiments on magma-filled cracks: Competition between external stresses and internal pressure. Earth Planet. Sci. Lett. 54, 1, 247–1, 261.

Wauthier, C., Cayol, V., Kervyn, F., d'Oreye, N., 2012. Magma sources involved in the 2002 Nyiragongo eruption, as inferred from an InSAR analysis. J. Geophys. Res.: Solid Earth 117, http://dx.doi.org/10.1029/2011JB008257.

Wright, T.J., Sigmundsson, F., Pagli, C., Belachew, M., Hamling, I.J., Brandsdóttir, B., Keir, D., Pedersen, R., Ayele, A., Ebinger, C., et al., 2012. Geophysical constraints on the dynamics of spreading centres from rifting episodes on land. Nat. Geosci. 5, 242–250.

Zia, H., Lecampion, B., 2020. Pyfrac: A planar 3d hydraulic fracture simulator. Comput. Phys. Comm. 255, 107368. http://dx.doi.org/10.1016/j.cpc.2020.107368.