# Nyamulagira's magma plumbing system inferred from 15 years of InSAR

Christelle Wauthier, Valérie Cayol, Mickael Poland, François Kervyn, Nicolas d'Oreye, Andrew Hooper, Sergey Samsonov, Kristry Tiampo, Benoit Smets

► To cite this version:

Christelle Wauthier, Valérie Cayol, Mickael Poland, François Kervyn, Nicolas d'Oreye, et al. Nyamulagira's magma plumbing system inferred from 15 years of InSAR. The Geological Society, London, Special Publications, 2013, 380 (1), pp.39-65. 10.1144/SP380.9. hal-03049852

HAL Id: hal-03049852

https://uca.hal.science/hal-03049852

Submitted on 10 Dec 2020

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
Nyamulagira's magma plumbing system inferred
from 15 years of InSAR

C. WAUTHIER¹,²,³*, V. CAYOL⁴, M. POLAND⁵, F. KERVYN², N. D'OREYE⁶,⁷,
A. HOOPER⁸,⁹, S. SAMSONOV¹⁰, K. TIAMPO¹¹ & B. SMETS²,⁷,¹²

¹GeMMe Unit, Georesources and Geo-Imaging, University of Liege, Liege, Belgium

²Earth Sciences Department, Royal Museum for Central Africa, Tervuren, Belgium

³Present address: Department of Terrestrial Magnetism, Carnegie Institution
of Washington, Washington, DC, USA

⁴Lab. *Magmas et Volcans, Université Blaise Pascal, UMR 6524, OPGC, Clermont-Ferrand, France*

⁵US Geological Survey, Hawaiian Volcano Observatory, Hawaii National Park, Hawaii, USA

⁶Departments of Geophysics/Astrophysics, National Museum of Natural History,
Walferdange, Luxembourg

⁷European Center for Geodynamics and Seismology, Walferdange, Luxembourg

⁸Department of Geoscience and Remote Sensing, Delft University of Technology,
Delft, The Netherlands

⁹Present address: School of Earth and Environment, University of Leeds, Leeds, UK

¹⁰Department of Earth Sciences, Western University, Ontario, Canada

¹¹Canada Centre for Remote Sensing, Ottawa, Canada

¹²Department of Geography, Vrije Universiteit Brussel, Brussels, Belgium

*Corresponding author (e-mail: cwauthier@dtm.ciw.edu)

**Abstract:** Nyamulagira, located in the east of the Democratic Republic of Congo on the western branch of the East African rift, is Africa's most active volcano, with an average of one eruption every 3 years since 1938. Owing to the socio-economical context of that region, the volcano lacks ground-based geodetic measurements but has been monitored by interferometric synthetic aperture radar (InSAR) since 1996. A combination of 3D Mixed Boundary Element Method and inverse modelling, taking into account topography and source interactions, is used to interpret InSAR ground displacements associated with eruptive activity in 1996, 2002, 2004, 2006 and 2010. These eruptions can be fitted by models incorporating dyke intrusions, and some (namely the 2006 and 2010 eruptions) require a magma reservoir beneath the summit caldera. We investigate inter-eruptive deformation with a multi-temporal InSAR approach. We propose the following magma plumbing system at Nyamulagira by integrating numerical deformation models with other available data: a deep reservoir (c. 25 km depth) feeds a shallower reservoir (c. 4 km depth); proximal eruptions are fed from the shallow reservoir through dykes while distal eruptions can be fed directly from the deep reservoir. A dyke-like conduit is also present beneath the upper southeastern flank of Nyamulagira.

The study and monitoring of active volcanic areas in Africa using ground-based methods can be jeop-
ardized by political tensions, security problems,
lack of sustainable infrastructures and difficult
accessibility. Until recently, the Lake Kivu area in the East African rift has suffered from a lack of ground-based geodetic measurements despite being frequently affected by earthquakes and volcanic eruptions. Indeed, two active volcanoes threaten

the population of the Kivu area: Nyamulagira and
Nyiragongo. As a result, remote-sensing tech-
niques, particularly interferometric synthetic aper-
ture radar (InSAR), provide useful and robust
observational tools for volcano monitoring and
hazard assessment in this area (d'Oreye et al. 2008).
Since 1996, SAR satellites have acquired images
of the Kivu area on a regular basis. InSAR provides
ground displacement maps, while multi-temporal

From: PYLE, D. M., MATHER, T. A. & BIGGS, J. (eds) *Remote Sensing of Volcanoes and Volcanic Processes: Integrating Observation and Modelling*. Geological Society, London, Special Publications, **380**.
http://dx.doi.org/10.1144/SP380.9 © The Geological Society of London 2013. Publishing disclaimer:
www.geolsoc.org.uk/pub_ethics
InSAR (MT-InSAR) can be used to overcome some of the limitations of conventional InSAR to detect potential inter-eruptive deformation and provide a time-series of surface displacements.

Nyamulagira is Africa's most active volcano, with an average of one eruption every 3 years since 1938 (Burt et al. 1994; Smets et al. 2010a), but its magma plumbing system is poorly known. An aseismic area lying between 3.5 and 7 km beneath the volcano is interpreted by Hamaguchi & Zana (1983) as a magma reservoir. This depth is consistent with tilt (Kasahara 1983), magnetotelluric sounding (Mishina 1983) and petrologic (Chakrabati et al. 2009; Head et al. 2011) studies. Wadge & Burt (2011) infer that the distribution of dykes feeding historic eruptions is indicative of a shallow reservoir (<10 km depth) that is larger than the caldera (c. 2 km radius) and slightly elongated in a NNW-SSE direction. Lava flows emitted from distal vents in 1904 and 1912, however, have petrologically primitive compositions (Denaeyer 1972; Aoki et al. 1985) that indicate a deeper (c. 20-30 km; Pouclet 1976) source and no shallow crystallization, suggesting that magma may bypass any shallow storage reservoirs on the way to the surface.

Modelling of geodetic data is a powerful tool to characterize the location and geometry of magma reservoirs and transport pathways. The aim of this paper was to study and model Nyamulagira InSAR data using the 3D mixed boundary elements method to better describe the magma plumbing system beneath the volcano. First, we present the geological setting of the Kivu area and the main characteristics of Nyamulagira eruptions. We then describe the InSAR data and processing methods, as well as our modelling procedure (including a method that provides a realistic distribution of overpressure, takes topography into account and is combined with non-linear inversions). Next, the combined modelling and inversion procedures are applied to InSAR observations of Nyamulagira eruptions in 1996, 2002, 2004, 2006 and 2010. We then investigate Nyamulagira's inter-eruptive behaviour using InSAR time-series analysis. Finally, we define the geometry and location of magma sources by integrating model and InSAR time-series results with other geological and geophysical data.

## Geological setting

Nyamulagira (3058 m a.s.l.) is a basaltic shield volcano located in the eastern part of the Democratic Republic of Congo (DRC), about 30 km NW of the city of Goma and 13 km NW of neighbouring Nyiragongo volcano (Fig. 1). Both volcanoes are part of the Virunga Volcanic Province

(VVP) in the western branch of the East African rift (EAR), which is characterized by an en-echelon succession of half-graben basins that are often filled by water to form lakes (Ebinger 1989a, b; Ebinger & Furman 2002/2003).

Cone and eruptive fissure distributions at Nyamulagira mainly follow a NNW-SSE-trending fracture network (Verhoogen 1948; Pouclet 1976; Smets et al. 2010a), as well as a NE-SW trend (Wadge & Burt 2011). The main NNW-SSE-trend may be due to the stress field generated by the interplay between the magmatic systems of Nyamulagira and Nyiragongo (Gudmundsson & Andrew 2007) or the reactivation of faults in the Precambrian basement (Smets et al. 2010b; Toombs & Wadge 2012).

Nyamulagira's summit is characterized by a c. 2 km-radius caldera that has walls up to 100 ms high and is crossed by numerous fractures. In the NE part of the caldera is a c. 400 m-wide pit crater (Fig. 2) that is about 90 m deep. A large depression, about 60 m deep, owing to a collapse during the 1938-1940 eruption is present in the SW half of the caldera.

The NNW-SSE-trending fracture network (Fig. 1) linking Nyamulagira and Nyiragongo crosses the caldera (Fig. 3a). On the upper part of the volcano's southeastern flank (Fig. 3b), the fracture network defines a graben which could be the surface expression of a dyke intrusion (Rubin & Pollard 1988; Rubin 1992), as observed along recently emplaced dykes in the EAR (Calais et al. 2008; Grandin et al. 2009).

Between 1958 and 1994 Nyamulagira activity was dominated by eruptions on the lower flanks, but eruptions after 1994 (in 1996, 1998, 2000, 2001, 2002, 2004, 2006 and 2010; Fig. 1) were mainly concentrated within the caldera and along the NNW-SSE fracture network on the upper flanks (Smets et al. 2010b). The 2011-2012 eruption of Nyamulagira occurred from a vent located c. 12 km east of the caldera (Fig. 1), close to that of the 1989 eruption. The eruption started on 6 November 2011, sent lava flowing northward and lasted until April 2012. Preliminary estimates made by the GORISK team (see GORISK Web site:http://www.ecgs.lu/gorisk/nyamulagira/current-activity/november-2011-eruption, 2012) suggest a deformed area of more than 250 km² and an erupted volume of at least 81.5 × 10⁶ m³.

Wadge & Burt (2011) identified two main types of eruptions based on their duration and location. The first class includes distal eruptions mainly before 1980 (Wadge & Burt 2011) or 1994 (Smets et al. 2010b) that occur NE and SW of the caldera and have long durations (>80 days) with no lava erupted in the caldera. The 2011-2012 eruption also belongs to this class of distal eruptions and might indicate the start of a new cycle of activity
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Geological setting of the Virunga Volcanic Province (VVP), western branch of the EAR, north of Lake Kivu. The rift extension direction in this area is about N110°E, with an estimated extension rate of 2.8 mm a⁻¹ (Stamps *et al.* 2008). The Goma–Gisenyi urban area is outlined in white. A dashed line highlights Nyamulagira's caldera rim. Nyamulagira lava flows erupted between 1996 and 2010 and are labelled. The white rectangle gives the extent of Figure 20.</figcaption></figure>

at Nyamulagira. The second class includes proximal eruptions that generate smaller volumes than distal ones (by about half), are of shorter duration (<80 days) and often occur in the caldera (especially since 1980). This class of eruption is approximately equally distributed between the upper northern and southern flanks of the volcano (Wadge & Burt 2011). The different eruption classes may be controlled by tectonic factors. Wadge & Burt (2011) suggested that eruptions fed by EAR-perpendicular dykes (oriented NW–SE) would need to have larger magma overpressures to produce the same output volumes as eruptions fed by EAR-parallel (NE–SW-oriented) dykes.

VVP lavas are basic, silica-undersaturated and K-rich (Kampunzu *et al.* 1982; Aoki & Yoshida

1983; Aoki *et al.* 1985). Such a composition results in very low-viscosity lavas that are able to flow for tens of kilometres at high velocities (Tazieff 1977). The Nyamulagira lava field covers over 1100 km² (Fig. 1) and contains more than 100 flank cones. Nyamulagira eruptions therefore pose a significant hazard to the surrounding environment and sometimes threaten the local populations.

## Methods

### Radar interferometry (InSAR)

InSAR is widely used to study deformation associated with a broad range of volcanic and seismic events (e.g. Burgmann *et al.* 2000; Lu *et al.* 2005;
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Nyamulagira's summit (viewed from the North) is characterized by a c. 2 km-radius caldera, outlined in white, which has walls up to 100 m high. The pit crater in the NE and depression in the SW part of the caldera are outlined in black (picture taken on 17 January 2010).</figcaption></figure>

Wright *et al.* 2006; Yun *et al.* 2006; Simons & Rosen 2007; Fournier *et al.* 2010; Hooper *et al.* 2012 and references therein). The method provides high-spatial-resolution, low-cost measurements of ground displacements with centimetre to sub-centimetre accuracy (Massonnet & Feigl 1998) over distances on the order of 100 km. Such data are valuable in areas where access is difficult or instrument cost is prohibitive, as is the case in the DRC. Since 1996, Nyamulagira has been

monitored by InSAR with data from the JERS, ERS-1/2, ENVISAT, RADARSAT-1/2, ALOS and TerraSAR-X satellites each being acquired for some part of the intervening years.

Conventional InSAR techniques suffer from several limitations. Dense vegetation covering most of the Kivu area causes strong decorrelation of the InSAR signal. This problem can be mitigated by using C-band ($\lambda \sim 5.6$ cm) interferograms having small temporal and perpendicular baselines (Fig. 4; d'Oreye *et al.* 2008), or L-band ($\lambda \sim 23$ cm) interferograms that can 'see' through some vegetation by employing a longer wavelength (Wauthier *et al.* 2009). Some areas remain incoherent, however, such as the vegetated slopes of the volcanoes and areas outside the rift valley (Fig. 4). Another limitation of conventional InSAR is contamination of the deformation signal by atmospheric delays or errors in the digital elevation model used to remove the topographic signal. These effects can be identified by comparing interferograms spanning different time intervals using a pair-wise logic method (Massonnet & Feigl 1995). The orbital errors are approximately planar (Hanssen 2001) and the interferograms shown in this study have thus been detrended to account for them. Finally, the estimation of the integer ambiguities in the phase known as 'phase unwrapping' can introduce 'phase jumps' that propagate over the entire scene. As we use the unwrapped data in our inversion procedure (see the section 'Non-linear inversion'), the unwrapped interferograms must be carefully checked for errors. However, new promising approaches directly invert the wrapped data by estimating the integer ambiguities simultaneously with the model parameters (Feigl & Thurber 2009; Hooper 2010; Ali & Feigl 2012).

Nyamulagira eruptions since 1996 are generally associated with localized ground displacements covering c. 50 km² on the upper flanks of the volcano (Colclough 2005; Wauthier *et al.* 2009; Wauthier 2011). An exception is the 2006 eruption, which occurred low on the southern flank of the volcano with no summit activity (Tedesco *et al.* 2007) and was associated with ground displacements over an area of c. 400 km² extending SE toward Nyiragongo (d'Oreye *et al.* 2008; Cayol *et al.* 2010; Cayol *et al.* in preparation).

An MT-InSAR approach can aid with investigations of inter-eruptive deformation. Indeed, MT-InSAR, which involves the processing of multiple SAR acquisitions over the same area, addresses some of the limitations of conventional InSAR. MT-InSAR algorithms can be grouped into two broad categories: persistent scatterer (PS) and small baseline (Hooper 2008). Each method is designed for a specific radar scattering mechanism. In the PS method, a pixel is dominated by a single

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** (a) Upper SE flank of Nyamulagira volcano (pictures were taken on 25 January 2010) showing a graben feature. Summit caldera rim is outlined in white. (b) Close-up of the graben feature.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Coherence images from ENVISAT ascending interferograms of the Nyamulagira area. The time spans ($B_t$) for the a-c coherence images are 35, 315 and 35 days, respectively. The perpendicular baselines ($B_p$) for the A, B and C coherence images are 77, 75 and 826 m, respectively.</figcaption></figure>

scatterer while the contributions from other scatterers are negligible. In contrast, a small baseline pixel includes distributed scatterers which have a phase with little decorrelation over short time periods. Here, we apply the 'StaMPS' ('Stanford Method for Permanent Scatterers') technique (Hooper 2008), which incorporates both the PS and small baseline approaches, to ENVISAT tracks that contain at least 20 images acquired between 2003 and 2010 over Nyamulagira. In the VVP, the InSAR time-series provide deformation maps with sub-centimetre accuracy even in areas that seem usually noisy in conventional C-band interferograms (Fig. 4).

### *Three-dimensional mixed boundary elements method*

We model surface deformation using the mixed boundary element method (MBEM; Cayol &

Cornet 1997), which takes into account realistic topographies and mechanical interactions between sources. The volcanic edifice is assumed to be linearly elastic, homogeneous and isotropic. Young's modulus and Poisson's ratio are fixed to values of 5 GPa and 0.25 respectively – values estimated from a seismic velocity profile of the VVP (Mavonga 2010) with Young's modulus corrected by a factor of 0.25 to account for the difference between static and dynamic modulus (Cheng & Johnston 1981; Wauthier et al. 2012).

When dykes reach the ground surface, they are assumed to be roughly quadrangles with their top part corresponding to the eruptive fissure and their bottom part corresponding to a straight line defined by six parameters (Fig. 5a) chosen to restrict the model search to mechanically plausible geometries (Fukushima et al. 2005, 2010). Three parameters–*Dip* (dip angle), *Shear* (angle between the dip direction and the line connecting the
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Modelling parameters. Three types of source model use the following geometrical parameters: a dyke that reaches the surface (a); a dyke that does not reach the surface (b); and an elliptical reservoir (c).</figcaption></figure>

midpoints of the dyke top and bottom) and Botelv (elevation of the bottom line midpoint with reference to sea-level) – define the location of the bottom line midpoint. Three other parameters – Botlen (length of the bottom of the dyke scaled to that of the top), Twist (horizontal angle between the top and bottom) and Botang (vertical angle of the bottom side) – give the position of the two end points of the dyke bottom. The top line location of the main quadrangular part of the dyke can also be inverted if needed by adding up to four extra parameters (northern and southern end coordinates in Fig. 6a). If the dyke is connected to several fissures at the surface, the along-dip length of the segmented parts, D_top, is added (Fig. 6b).

Dykes not connected to the ground surface are assumed to have a simpler quadrangular shape, with a straight top and bottom (Fig. 5b). The top, assumed to be horizontal, is defined by three parameters, Toplength (length), D_top (average depth below the ground surface) and strike (angle with respect to the north direction). The bottom is defined by four additional parameters: Dip (dip angle), Botelv (elevation of the bottom midpoint

above sea-level), Botlen (length of the dyke bottom scaled to that of the top) and Botang (the vertical angle of the dyke bottom).

Magma reservoirs are represented by ellipsoidal sources defined by eight geometrical parameters (Fig. 5c). The parameters midx, midy and midz define the position of the ellipse centre, S1 and S2 define the lengths of the two axes (with S1 > S2), Rotplane gives the angle of rotation for the ellipse in the xy plane (counter-clockwise from the x-axis) and Dipdir is the dip direction (counter-clockwise from the x-axis) of the dip angle Dipₑ (measured from horizontal; a horizontal ellipse will thus have a Dipₑ = 0; Fukushima 2006).

In contrast to simple analytical solutions for sources in an elastic half-space, for example Okada (Okada 1985) for rectangular dislocations and Mogi (Mogi 1958) or Yang (Yang *et al.* 1988) for spherical or ellipsoidal sources, boundary conditions are imposed on each type of meshed source defined above. The prescribed boundary conditions are tractions that are assumed to be perturbations of an initial state of stress. Traction are zero at the ground surface and are equal to overpressures in

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** (a) The top of the main quadrangular part of the model dyke is parallel to the curve going through the fissure midpoints, and the upper surface of the dyke is divided into segments that intersect each eruptive fissure. (b) Schematic figure explaining the parameterization of the segmented part of a model dyke. D_top is the along-dip depth from the surface of the segmented parts. Modified after Fukushima (2006).</figcaption></figure>
dykes or reservoirs (Cayol & Cornet 1997). With the goal to constrain the stress field in the lower duc- tile Icelandic crust, Hooper et al. (2011) include a lithostatic pressure gradient to model a deep dyke. They defined this gradient by assuming a density difference of 250 kg m⁻³ between the magma and the surrounding rock. A study by Fukushima et al. (2010), however, showed that pressure gradients cannot be constrained in the upper crust at Piton de la Fournaise volcano (Reunion Island). At Ny- mulagira, preliminary models (Cayol et al. 2010; Wauthier 2011) showed that the dykes are shallow. Furthermore, we do not have an accurate a priori knowledge of the stress field (the overpressure) of the shallowest part of the edifice. In this study, con- stant dyke overpressures are therefore assumed. This boundary condition leads to a non-uniform dis- tribution of displacements. The source geometries are inverted simultaneously with the tractions, lead- ing to a minimum of seven parameters per source.

Topography is meshed by planar triangular elements. In order for the limited extension of the mesh to have a negligible influence on computed displacements (<2%), the size of the topographic mesh is 10 times larger than the maximum horizon- tal dimension of the deformation source (Cayol 1996). The density of the topographic mesh is fine close to the deformation sources, where displace- ment gradients are the largest, and becomes coarser with distance. A previous study by Wauthier (2011) showed that the node density close to the eruptive fissure(s) is the most critical parameter for model accuracy. For the Nyamulagira 1996 erup- tion, for example, the mesh (Fig. 7a) has a radius of 30 km, uses an interval of 50 m between the nodes adjoining the eruptive fissures, and is com- posed of 2092 elements and 1048 nodes overall.

Source boundaries also are meshed with planar triangular elements. Cayol & Cornet (1997) showed that coarse meshing of sources induces an overestimation of the amplitude of surface defor- mation, leading to an underestimation of source overpressure (Fukushima et al. 2005). A good com- promise between calculation time and model accu- racy must be defined on a case-by-case basis (Fukushima et al. 2005; Wauthier 2011). As an example, the node density for the source mesh is equal to a 200 m interval between the nodes for the 1996 Nyamulagira eruption.

Non-linear inversion

Following Fukushima et al. (2005), we simul- taneously invert for the dyke geometry, location and uniform overpressure. As surface displace- ments are a non-linear function of the source par- ameters, a non-linear inversion, based on a nearest neighbourhood inversion algorithm (NA), is used. The algorithm consists of two stages: search (Sambridge 1999a) and appraisal (Sambridge 1999b). A misfit function quantifies the discrep- ancy between observed and 3D-MBEM modelled displacements (Fukushima et al. 2005):

$$
\chi^2 = (u_o - u_c)^T C_d^{-1} (u_o - u_c) \quad (1)
$$

where C_d is the data covariance matrix and u_o and u_c, formed by concatenating data from the differ- ent viewing geometries, are vectors of observed and modelled line-of-sight (LOS) displacements, respectively.

To make the misfit computation numeri- cally manageable, u_o is constructed by subsampling the unwrapped LOS displacements (Fig. 7b).

<figure><img src="image_8.png" /><figcaption>Fig. 7. (a) Topographic mesh overlain on the shaded relief map (based on the SRTM DEM) used for modelling deformation from InSAR for the 1996 eruption. (b) Unwrapped JERS Interferogram (19 October 1996 to 28 February 1997) spanning the 1996 Nyamulagira eruption. (c) Data points subsampled from the JERS InSAR data set covering the 1996 eruption. Nyamulagira's summit caldera is outlined in white.</figcaption></figure>
Fukushima et al. (2005) showed that, for dykes at Piton de la Fournaise volcano, the subsampling method does not affect the result; hence a simple circular-grid subsampling method was chosen. Data, excluding pixels with a coherence of less than 0.15, are subsampled in such a way that the number of subsampled points is denser close to the deformation sources and coarser further away. After subsampling, the 5.5 million points of an original JERS-1 descending interferogram covering the 1996 Nyamulagira eruption were reduced to 1966 points (Fig. 7c). The offset of each unwrapped dataset is determined for each forward model calculated during a nearest-neighbour inversion by minimizing the misfit following Wauthier et al. (2012). Note that, in the present study, we show the observed, modelled and residual interferograms as the wrapped phase at the original sampling, to better depict the detail of the deformation field.

InSAR data contain spatially correlated phase noise caused by orbit errors as well as propagation of signal through the atmosphere. Statistical studies of atmospheric noise are based on computations of either the autocorrelation or the covariance function (Tarantola 1987). Both can be fitted by exponential functions of the form:

$$C(r) = \sigma_d^2 \exp\left(\frac{-r}{a}\right) \quad (2)$$

where r is the distance between two pixels, σ (2/d) the variance of the noise, and a is the correlation length (Fukushima et al. 2005; Parsons et al. 2006). Typical values determined from the residuals of preliminary inversions are 10⁻⁴ and 300 m for the variance and correlation distance, respectively, and are used in our initial inversions. These values are modified if they are too far from the values determined from the residuals of the best-fit model and the inversions are re-run (Fukushima et al. 2010). Indeed, Fukushima et al. (2005) showed that taking the data correlation into account has a significant effect on the best-fitting model, but that doubling the correlation length has a negligible effect on the results.

In the search stage of the inversion, the model space is explored within predefined parameter bounds. First, n1 (here 200) initial acceptable models are chosen randomly and the corresponding misfits are calculated. In the next iteration, n new models are generated in the neighbourhood (Voronoi cell) of the n lowest-misfit models, and misfits for these new models are calculated. Iterations continue until the misfit is no longer lowered significantly. Increasing n will make the search slower but will allow a larger parameter space to be explored (Sambridge 1999a). Fukushima et al. (2005) showed that values of 10–70 for

n lead to robust results. Since the inversions of Nyamulagira InSAR data are computationally extensive, n is taken to be 10 or 30.

In the appraisal stage (Sambridge 1999b), posterior probability density functions are calculated using misfit values computed during the search stage. Model statistical characteristics, such as the mean model and marginal posterior probability density functions (PPDs), are estimated. Confidence intervals for the model parameters can be obtained from the one-dimensional marginal PPDs, while trade-offs between pairs of model parameters are given by the two-dimensional marginal PPDs. The applicability of the method has been confirmed by synthetic tests (Fukushima et al. 2005).

To compare inversion results with a physically meaningful characteristic, we introduce an RMS error, expressed in centimetres, defined as:

$$RMS = 100 * \left( \sqrt{\frac{(\mathbf{u}_o - \mathbf{u}_c)^T (\mathbf{u}_o - \mathbf{u}_c)}{\text{npts}}} \right) \quad (3)$$

where npts is the number of subsampled points.

## Case study: the December 1996 Nyamulagira eruption

On 1 December 1996, a flank eruption began along two NNW-oriented fissures located 2.5 km NW of the caldera rim (Fig. 8). The eruption lasted at least 4 days (BGVN 1996a) and emitted an estimated lava volume of 49.8 ± 16.6 × 10⁶ m³ (Smets et al. 2010a). The eruption was preceded in September 1995 by an intense period of seismicity (BGVN 1996b), as well as several hours of low amplitude tremor following a series of tectonic earthquakes between 15 November and 2 December 1995. An increase in seismicity might have occurred again in August 1996, four months before the eruption. However, we would like to point out that the BGVN reports referenced in this study are incomplete and that the seismic data are poorly constrained as installing and maintaining ground-based networks in the Kivu area is a challenging task.

An L-band JERS descending interferogram covers the period from 19 October 1996 to 28 February 1997 (Fig. 8). This interferogram spans the 1996 eruption but may also include both pre-eruptive (between 19 October and 1 December 1996) and post-eruptive (between 17 December 1996 and 28 February 1997) deformation. The interferogram shows a c. 50 km² extension of the deformation signal associated with the eruption on the volcano’s NW flank. In Figure 8, signal A, east of the eruptive fissures, corresponds to approximately
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Wrapped JERS interferogram from Figure 7b with the 1996 lava flows and the eruptive fissures (two bold lines) shown. Signals a–c are discussed in the text. One colour cycle represents 11.8 cm of line of sight (LOS) range change. Height of ambiguity ($H_{amb}$, the altitude difference that generates one topographic-related fringe) is 380 m. Nyamulagira's summit caldera is outlined with a dashed white line.</figcaption></figure>

three fringes of range decrease (35.4 cm of displacement towards the satellite), while signals B and D correspond to small range increases of less than one fringe (<11.8 cm of displacement away from satellite). Finally, signal C corresponds to a small range decrease of less than one fringe.

The asymmetric pattern of deformation is consistent with dyke intrusions associated with the two eruptive fissures. Signal C could be induced by dyke intrusion or by an inflating magma reservoir beneath the SW depression within Nyamulagira's caldera, as suggested by Kitagawa et al. (2007) (note that the deflating magma reservoir suggested by Toombs & Wadge (2012) has to be ruled out as these authors misinterpreted signal C as subsidence because they incorrectly treated the JERS interferogram as being acquired in an ascending geometry; Geoff Wadge, pers. comm. August 2012). Kitagawa et al.'s (2007) preliminary forward model for the deformation pattern associated with the 1996 eruption includes a dyke 0.6 km high, 1.4 km long, reaching a depth of 1 km, and with an average opening of 1.5 m.

We used the combined 3D-MBEM and NA inversion procedure and found that the deformation can be explained with a single vertical dyke reaching the ground surface through two en-echelon eruptive fissures. The best-fitting dyke (Fig. 9) is

vertical and c. 3 km high. The two en-echelon fissures are connected 400 m beneath the surface to the main contiguous part of the dyke. The top of the main dyke is laterally more extensive than the surface eruptive fissures, but the dyke tapers to a smaller length at its base ($Botlen = 0.26$). The modelled dyke overpressure is 1.3 MPa, which is consistent with typical overpressures found for other Nyamulagira and Nyiragongo eruptions (Wauthier 2011; Wauthier et al. 2012).

The bottom of the best-fitting dyke is located 3.4 km beneath the summit of the volcano, which is consistent with the presumed depth of Nyamulagira's magma chamber (from 3.5 to 7 km depth). We tested such an inflating reservoir below the dyke in our inversions, but the result was zero overpressure for the reservoir, indicating that any pressure change in a shallow reservoir pressure was too small to be detected. We therefore cannot constrain the location of any magma reservoir beneath Nyamulagira from these data.

The model dyke has a volume of $5.56 \times 10^6 \text{ m}^3$ and induces an average opening of 0.75 m with a maximum surface opening of 1 m, which is consistent with typical values of Nyamulagira eruptive fissures (Pouclet 1976). InSAR data are not coherent over the entirety of Nyamulagira's flanks, so the dyke could be larger than proposed here. The dyke dimensions are about 3.4 km high, 2.2 km long and with 0.75 m opening. This is a larger extension and a smaller opening than the preliminary forward model from Kitagawa et al. (2007).

## Models of other Nyamulagira eruptions

We applied the same combined 3D-MBEM and NA inversion procedures to Nyamulagira eruptions captured by InSAR in 2002, 2004, 2006 and 2010. Unfortunately, InSAR data covering the 1998, 2000 and 2001 eruptions are too noisy to allow for modelling. As established by Wauthier (2011) and Toombs & Wadge (2012), most of the recorded co-eruptive InSAR deformation at Nyamulagira can be explained by dyke intrusions.

### *The 1998, 2000 and 2001 eruptions*

**The 1998 eruption.** On 17 October 1998, a flank eruption began and lasted for at least 8 days (BGVN 1998). The eruptive fissure, starting from the caldera rim, opened on the NW flank of the volcano and emitted an estimated lava volume of $69.2 \pm 23.1 \times 10^6 \text{ m}^3$ (Smets et al. 2010a). Two ascending interferograms, from the C-band RADARSAT-1 and ERS (Toombs & Wadge 2012) satellites, contain some deformation associated with the eruption. Owing to the large time
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** (a) Observed (left), modelled (middle) and residual (right) interferograms based on a model of a single dyke intrusion associated with the 1996 eruptive fissures (two bold lines). Nyamulagira's summit caldera is outlined with a dashed white line. (b) Geometry of the best-fit single dyke model. The view with contour lines is a map view. The dashed lines in the map view indicate the locations of the vertical cross-sections shown to the right of, and below, this view. The view to the right is the north-vertical cross-section and the view below is the east-vertical cross-section.</figcaption></figure>

spans, however, the coherence is low (Fig. 10a). The asymmetric pattern of deformation is consistent with one or more dyke intrusions, but the noise is too high and the signal too unclear for in-depth modelling analysis.

*The 2000 eruption.* A flank eruption started on 27 January 2000 and lasted approximately 14 days (BGVN 2000). The eruptive fissure extended from the caldera rim towards the southeastern flank over a distance of about 1.3 km and effused an estimated lava volume of $47.3 \pm 15.8 \times 10^6 \text{ m}^3$ (Smets et al. 2010a).

A RADARSAT-1 descending interferogram shows some deformation associated with the 2000

eruption (Fig. 10b), but coherence is low (similar to the 1998 eruption interferograms). A range increase, centred on the pit crater inside Nyamulagira's caldera, is the only visible ground displacement. This signal could be associated with reservoir deflation or dyke intrusion. Owing to lack of other coherent data and the large time span, however, we cannot rule out that this signal is contaminated by compaction of older lava flows or by atmospheric artefacts.

*The 2001 eruption.* This long-lived flank eruption started on 6 February and lasted 54–62 days (BGVN 2001a, b, c). The lava flowed both north and south from the caldera rim, with the formation
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** (a) Ascending RADARSAT-1 wrapped interferogram (19 February 1997 to 5 November 1998; incidence angle $c$ 44°, $H_{amb} = 1480$ m) showing deformation associated with Nyamulagira's 1998 eruption. (b) Descending RADARSAT-1 wrapped interferogram (26 October 1998 to 25 January 2000, incidence angle $c$ 35°, $H_{amb} = 27$ m) showing deformation associated with Nyamulagira's 2000 eruption. (c) ERS-2 ascending wrapped interferogram (6 September 2000 to 3 July 2002, incidence angle $c$ 23°, $H_{amb} = 147$ m) showing deformation associated with Nyamulagira's 2001 eruption. In all images, lava flows are mapped and the eruptive fissures are shown as two bold lines. Nyamulagira's summit caldera is outlined with a dashed white line.</figcaption></figure>

of two new cones and an estimated lava volume of 140.1 ± 46.7 × 10⁶ m³ (Smets et al. 2010a).

We were not able to find a pair of SAR images suitable to capture the deformation induced by this eruption without including the major Nyiragongo 2002 dyking event (Wauthier et al. 2012). Furthermore, the Nyamulagira's co-eruptive signal (Fig. 10c) is of limited extent on the incoherent volcano flanks (Wauthier et al. 2012; Toombs & Wadge 2012); the deformation is therefore very difficult to interpret. Its limited extension is somewhat surprising given the long eruption duration and large emitted volume.

## The 2002 eruption

The 2002 eruption started on 25 July and lasted from 15 to 63 days (BGVN 2002a, b, 2003). Fissures opened in the caldera and extended north outside the caldera boundaries (Fig. 11), and the estimated lava flow volume was 56.7 ± 18.9 × 10⁶ m³ (Smets et al. 2010a). An increase in seismicity was observed at Nyamulagira from the end of February 2002 until the eruption, and a period of low seismic activity, characterized by shallow long-period (LP) earthquakes, persisted for about a year following the eruption's end (Mavonga et al. 2006).

The 2002 eruption was spanned by SAR data from ascending ERS-2 orbits and ascending and descending RADARSAT-1 orbits (Fig. 11). NW of the eruptive fissure, a signal (labelled A in Fig. 11) corresponding to a range decrease (motion

towards the satellite) is visible in all interferograms. On the RADARSAT-1 ascending interferogram, 15 narrow fringes corresponding to approximately 42 cm of range decrease are apparent (signal A). The sense of signal A switches to range increase (motion away from the satellite) at signal D, located north of the eruptive fissure.

Signal A is also visible on the RADARSAT-1 descending interferogram, with about three fringes of range increase. The sense of motion of the deformation changes to the SE (black arrow in Fig. 11c) and becomes three narrow fringes of range decrease located west of the eruptive fissure (signal C). An additional deformation signal (B) east of the eruptive fissure and SE of the caldera corresponds to a small range decrease (1.5 fringes).

The A and D signals in ascending interferograms and the B signal in the descending interferogram are probably associated with opening of the eruptive fissure connected to one or more dykes. After a check by generating coherent interferograms with other RARDARSAT descending pairs, we can rule out that signal C is an atmospheric effect, and the deformation could indicate another source, such as an inflating sill-like reservoir or another dyke which did not reach the surface.

Colclough (2005) previously identified the deformation signal linked with the 2002 eruption on ERS interferograms and concluded that a simple Mogi source was not able to reproduce the deformation pattern. Owing to the asymmetry of the pattern associated with the eruptive fissure, the
<figure><img src="image_12.png" /></figure>

**Fig. 11.** Three wrapped interferograms showing deformation associated with Nyamulagira's 2002 eruption. The 2002 lava flows are mapped and the eruptive fissure is shown as a bold line. See text for interpretations of signals labelled A–C. (a) ERS-2 ascending interferogram (3 July 2002 to 29 January 2003, incidence angle $c. 23^\circ$, $H_{amb} = 3900$ m). (b) Ascending RADARSAT-1 interferogram (17 February 2002 to 15 October 2002, incidence angle $c. 44^\circ$, $H_{amb} = 36$ m). (c) Descending RADARSAT-1 interferogram (3 March 2002 to 11 September 2002, incidence angle $c. 35^\circ$, $H_{amb} = 23$ m). Nyamulagira's summit caldera is outlined with a dashed white line.

most likely source is a dyke connected to the surface through the mapped eruptive fissure and trending north-south (Toombs & Wadge 2012).

The three interferograms shown in Figure 11 were inverted to model deformation associated with the 2002 eruption. The geometry of the best-fitting model is a dyke that is longer at depth than the eruptive fissure, dipping east at shallow levels, and dipping 40° west from 1 km beneath the volcano's summit to its base at 3.5 km depth (Fig. 12). The depth to the dyke bottom is consistent with the estimated depth of Nyamulagira's magma reservoir. The dyke volume is $27 \times 10^6$ m$^3$, which is about half of the estimated erupted volume (Table 1), and the modelled dyke overpressure is about 0.4 MPa.

## The 2004 eruption

In 2004, an eruption started on 8 May and lasted from 20 to 23 days on the volcano's NW flank, 1.9 km from the caldera rim (Fig. 13). The estimated lava flow volume was $68.5 \pm 22.8 \times 10^6$ m$^3$ (Smets et al. 2010a). At the beginning of the eruption, volcanic activity was also reported to have occurred inside the caldera (BGVN 2004a).

Seismic activity (including LP earthquakes attributed to the movement of magma in a deep conduit) increased around Nyamulagira about 10 months prior to the 2004 eruption (Mavonga et al. 2006). A swarm of LP earthquakes, accompanied by tectonic events possibly associated with the EAR, began about four months prior to the eruption.

**Fig. 12.** (a) Observed (left), modelled (middle) and residual (right) interferograms based on a model of a single dyke intrusion associated with the 2002 eruptive fissures (bold line). Nyamulagira's summit caldera is outlined with a dashed white line. (b) Geometry of the best-fit single dyke model. Views are as in Figure 9.
NYAMULAGIRA'S MAGMA PLUMBING SYSTEM

(a)
North (km)

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

(b)

East (km)

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

Elevation (m)

Elevation (m)

<figure><img src="image_15.png" /></figure>
<table><caption>**Table 1.** Characteristics of modelled Nyamulagira eruptions</caption><thead><tr><th rowspan="2">Fissure location</th><th rowspan="2"></th><th>1996</th><th>2002</th><th>2006*</th><th>2010</th></tr><tr><th>North flank</th><th>North flank and caldera</th><th>SE flank</th><th>SE flank and caldera</th></tr></thead><tbody><tr><td>(P<sub>0</sub>)<sub>dyke1</sub></td><td>(MPa)</td><td>1.3</td><td>0.4</td><td>0.2</td><td>1</td></tr><tr><td>(P<sub>0</sub>)<sub>dyke2</sub></td><td>(MPa)</td><td>–</td><td>–</td><td>–</td><td>3</td></tr><tr><td>Reservoir centre</td><td>(km)</td><td>–</td><td>–</td><td>744, 9843</td><td>744.7, 9843.4</td></tr><tr><td>Reservoir dimensions</td><td>(km)</td><td>–</td><td>–</td><td>1 × 1</td><td>1.7 × 1.2</td></tr><tr><td>Reservoir depth</td><td>(km below caldera)</td><td>–</td><td>–</td><td>1.5</td><td>3.7</td></tr><tr><td>(ΔV)<sub>reservoir</sub></td><td>10<sup>6</sup> m<sup>3</sup></td><td>–</td><td>–</td><td>+1.5</td><td>–2</td></tr><tr><td>(V)<sub>dyke1</sub></td><td>10<sup>6</sup> m<sup>3</sup></td><td>5.6</td><td>27</td><td>1.2</td><td>1.1</td></tr><tr><td>(V)<sub>dyke2</sub></td><td>10<sup>6</sup> m<sup>3</sup></td><td>–</td><td>–</td><td>–</td><td>2.7</td></tr><tr><td>(V)<sub>intruded</sub></td><td>10<sup>6</sup> m<sup>3</sup></td><td>5.6</td><td>27</td><td>2.7</td><td>3.8</td></tr><tr><td>(V)<sub>extruded</sub></td><td>10<sup>6</sup> m<sup>3</sup></td><td>49.8 ± 16.6</td><td>56.7 ± 18.9</td><td>44.2 ± 14.7</td><td>46</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Three wrapped interferograms showing deformation associated with the 2004 Nyamulagira eruption. Arrows show the LOS direction. The 2004 lava flows are mapped and the eruptive fissure is shown as the bold light line, while other Nyamulagira fissures are shown as dark bold lines. See the text for the interpretations of signals labelled A–C. (a) ENVISAT ascending interferogram (24 March 2004 to 2 June 2004, incidence angle c. 18°, H<sub>amb</sub> = 49 m). (b) Ascending RADARSAT-1 interferogram (26 March 2004 to 30 June 2004, incidence angle c. 44°, H<sub>amb</sub> = 26 m). (c) Descending RADARSAT-1 interferogram (16 March 2004 to 7 August 2004, incidence angle c. 35°, H<sub>amb</sub> = 30 m). Nyamulagira's summit caldera is outlined with a dashed white line.</figcaption></figure>

<footnote>*Denotes results from Cayol *et al.* (2010).</footnote>
About two months before the eruption, the number of deep LP events decreased and the number of shallow LP earthquakes increased. A major seismic swarm occurred on 4–6 May and a tectonic earthquake was felt in the VVP north of Lake Kivu about 4 h before the eruption started. About 2.5 h before the eruption, a high-amplitude, long-period earthquake occurred north of the VVP. After the eruption began, seismic activity was dominated by shallow (0–5 km) LP events (Mavonga et al. 2010).

Unusual fumarolic activity was observed inside Nyamulagira's caldera on 2 May (seen from Goma). Once the eruption started, a lava lake was present inside the pit crater within the caldera until at least 12 May. Four new cones (each 30–50 ms high) were built on the NW flank of the volcano during the eruption (BGVN 2004a, b, 2006).

The 2004 eruption was spanned by data from ENVISAT ascending orbits and RADARSAT-1 ascending and descending orbits (Fig. 13). North of the caldera is a small-gradient deformation signal (labelled A in Fig. 13), indicating a range increase in the ascending orbits and a range decrease in the descending orbits. Within the NW part of the caldera is a large deformation signal with a high gradient (signal B), and a similar signal best viewed in the ascending ENVISAT interferogram (Fig. 13a) is present just outside the NW margin of the caldera (signal C). SE of the caldera, a small deformation signal (D) indicates a range increase in the RADARSAT ascending interferogram (also present, but of small magnitude, on the descending interferogram as well).

Signal A may originate from the northern edge of the eruptive fissure, and thus be associated with an underlying dyke intrusion. The B and C signals probably correspond to a shallower phenomenon beneath the caldera, such as a dyke or sill intrusion. Signal D is probably an atmospheric effect, since it follows the volcano slopes and is not identified on other independent interferograms.

We modelled the deformation north of the caldera (signal A) with a shallow (bottom at 2.8 km below Nyamulagira's summit) intrusion (Fig. 14) dipping 11° to the west – almost a sill – with an overpressure of about 3 MPa. Deformation within and adjacent to the caldera (B and C signals) and on the southern flank (signal D), however, cannot be explained by such a model, and more data are required to constrain possible sources active during the 2004 eruption. The current dyke model is considered too preliminary and therefore will not be included in further discussion.

## *The 2006 eruption*

In 2006, an eruption started on 27 November with the opening of a 5 km-long eruptive fissure south

of the caldera rim, located half-way between Nyamulagira and Nyiragongo (Fig. 15). The lava flowed SW toward Lake Kivu, and the estimated lava volume was $44.2 \pm 14.7 \times 10^6 \text{ m}^3$ (Smets et al. 2010a). The eruption lasted between 8 and 23 days.

Several LP seismic swarms, with deep (20–30 km) and shallow (0–5 km) events separated by an aseismic zone, occurred beneath the volcano about a year prior the eruption. The aseismic zone is interpreted as hosting a magma reservoir (Mavonga et al. 2010), similar to the reservoir found by Hamaguchi & Zana (1983), while the LP events could be generated by flow of magma into Nyamulagira's shallow reservoir from a deep source. The shallow LP events were only recorded about one to two months before the eruption, probably because they were associated with the filling of a shallow reservoir (Mavonga et al. 2010). A magnitude 6.8 earthquake and numerous aftershocks (USGS NEIC catalogue) struck the Lake Tanganika area (located c. 500 km south of Nyamulagira volcano) in December 2005. Following this event, sustained seismic activity characterized by long-period earthquakes, interpreted as magma movement at depth, was recorded at Nyamulagira (BGVN 2007).

This 2006 eruption was spanned by ENVISAT interferograms from three different look angles, two from descending orbits and one from an ascending orbit (Fig. 16a). Only one interferogram from the ascending orbits had spatial and temporal baseline conditions that acted to minimize coherence loss SW of the volcano. This interferogram shows a signal over an area of 200 km² to the east, north and south of the eruptive fissure (Fig. 15). In contrast, interferograms from descending orbits only show displacements over a 25 km² area east and north of the fissures (Cayol et al. 2010).

In the ascending orbits, signal A in Figure 15, located SE of the caldera and NE of the eruptive fissures, corresponds to a range increase (motion away from the satellite), while signals B and C, located north and south of the eruptive fissures, respectively, correspond to range decreases (motion towards the satellite). Although they are separated by incoherent areas, signals B and C are consistent with one or more dyke intrusions associated with the eruptive fissures. The unusually large area of the deformation pattern suggests a large source dimension compared with the size of the eruptive fissures at the surface. This extensive deformation pattern differs from the ground displacements related to past eruptions captured by InSAR (described above). Comparison of the petrological characteristics of the 2006 lava to previously erupted lava flows on Nyamulagira's south flank
(a)

ENV. ASC.

RSAT ASC.

RSAT DESC.

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

DATA

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

MODEL

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

RESIDUAL

(b)

East (km)

<figure><img src="image_20.png" /></figure>
<figure><img src="image_21.png" /></figure>

**Fig. 15.** ENVISAT ascending interferogram (31 October 2006 to 5 December 2006, incidence angle $c$. 43°, $H_{amb} = 102$ m) spanning the 2006 Nyamulagira eruption, with the 2006 lava flows mapped and the eruptive fissures shown as a bold line. Nyamulagira's summit caldera is outlined with a dashed white line. See text for explanation of signals labelled A–C.

shows no significant compositional differences (Elisabet M. Head, pers. comm. August 2008).

Preliminary results (Cayol *et al.* 2010) demonstrate that a single dyke fits the deformation patterns poorly, and two sources are needed to explain the measured displacements. One of the sources corresponds to a dyke located between Nyamulagira and Nyiragongo, which is aligned with the axis linking the two volcanoes. The other source, which is responsible for displacements measured north of the eruption site, can be explained by an inflating sill-like source or a spherical reservoir (Cayol *et al.* 2010). The descending interferograms (Fig. 16), however, span a long time interval and could thus include pre- or post-eruptive deformation.

## *The 2010 eruption*

The 2010 eruption began on 2 January at 1:07 am (local time, UTC +2) according to seismic data

from the Goma Volcano Observatory monitoring network. Lava erupted from the pit crater located in the northeastern part of the caldera and along an ENE-WSW-trending fracture in the southern part of the caldera (Figs 17 & 18). A 600 m-long fissure also opened on the SSE flank of the edifice, about 1.5 km from the caldera rim, building a 150 m-high cone. The main flows from the flank fissure were 10.5-11 km long and partly covered portions of the 2006 lava flow, as well as the southern part of the 2001 lava flow (Fig. 1). The eruption lasted until 27 January, although the activity in the caldera probably stopped on 9 January (Smets *et al.* in preparation). The estimation volume for the 2010 eruption was $53 \pm 17.9 \times 10^6 \text{ m}^3$ (Smets *et al.* 2010a).

There was no noticeable increase in seismic activity in the months preceding the January 2010 eruption. According to D. Tedesco (pers. comm. January 2010), both tectonic earthquakes and tremor occurred simultaneously with the emission of lava at the surface. This observation suggests that magma was already close to the surface before the eruption. This, however, contrasts with the observed pre-eruptive deformation identified on time-series (see the section 'Complementary MT-InsAR results'). Ground displacements related to the January 2010 eruption were captured by the ENVISAT, RADARSAT-2, and ALOS satellites (Fig. 17 and Table 2).

In Figure 18, signal A, east of the northern eruptive fissure, corresponds to a range increase in the descending interferograms, but is not clearly visible in the ascending interferograms. Signals B and C, SE of the northern fissure and east of the southern fissure, respectively, correspond to range decreases in descending interferograms and range increases in ascending interferograms. Finally, signal D, which is only visible in descending interferograms, corresponds to a small range increase.

Signals A–C are consistent with one or more dyke intrusions associated with the eruptive fissures, as the sense of ground displacement changes depending on whether the interferograms are ascending or descending. This change in the sense of displacements is characteristic of horizontal deformation, such as that associated with dyke emplacement. Signal D could be caused by interactions between sources, another dyke that did not reach the surface or fault movement. Because the 2010 eruption was captured by numerous independent interferograms, we can constrain most ground deformations to have occurred between 29

**Fig. 14.** (a) Observed (left), modelled (middle) and residual (right) interferograms based on a model of a single dyke intrusion associated with the 2004 eruptive fissures (bold line). Nyamulagira's summit caldera is outlined with a dashed white line. (b) Geometry of the preliminary best-fit single dyke model. Views are as in Figure 9.
<figure><img src="image_22.png" /><figcaption>**Fig. 16.** (a) Observed (left), modelled (middle) and residual (right) interferograms based on a model of a dyke intrusion and an inflating reservoir associated with the 2006 eruptive fissures (bold line). Nyamulagira's summit caldera is outlined with a dashed white line. (b) Geometry of the preliminary model, adapted from Cayol *et al.* (2010). Views are as in Figure 9.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Fig. 17.** Wrapped interferograms covering the 2010 Nyamulagira eruption. Note that i#, F21N and UF denote ENVISAT beam number, RADARSAT-2 Fine mode and beam number and RADARSAT-2 Ultra-Fine mode, respectively. Eruptive fissures are shown as bold lines and Nyamulagira's summit caldera is outlined with a dashed white line. See Table 1 for data description.</figcaption></figure>

December 2009 and 10 January 2010, the latter date corresponding to the end of volcanic activity inside the caldera.

Eight interferograms, including six from ENVISAT (two ascending and four descending), one from RADARSAT-2 (descending) and one from ALOS (ascending), were inverted to model the 2010 eruption. The best-fitting model (Fig. 19) includes two dykes with two independent overpressures and a deflating sill-like reservoir below Nyamulagira caldera. The northern dyke dips 80°

eastward, has a volume of $1.1 \times 10^6 \text{ m}^3$, and reaches a depth of 2.6 km beneath the summit of the volcano, while the southern dyke dips 70° westward, has a volume of $2.7 \times 10^6 \text{ m}^3$, and reaches a depth of 3.7 km. Both bottom depths are consistent with the estimated depth of the top of the Nyamulagira magma reservoir. The northern dyke overpressure is three times lower (1 MPa) than that of the southern dyke (3 MPa). The RMS error for the best-fitting two dykes and deflating sill model is 1.8 cm, suggesting reasonable model residuals (Fig. 19a).
<figure><img src="image_24.png" /><figcaption>**Fig. 18.** (a) ENVISAT descending wrapped interferogram spanning 10 December 2009 to 14 January 2010, i7 beam mode (incidence angle c. 38°), and (b) ascending wrapped interferogram spanning 29 December 2009 to 2 February 2010, i7 beam mode (incidence angle c. 37°); both images show deformation patterns associated with the main phase of the 2010 eruption. Eruptive fissures are shown as bold lines and Nyamulagira's summit caldera is outlined with a dashed white line. See text for explanation of signals labelled A–D.</figcaption></figure>

Because of the two distinct dyke directions and overpressures, the two dykes probably reflect different emplacement mechanisms under different stress conditions governing the caldera and the SE flank. The northern dyke bottom is about 1 km shallower than the modelled sill, which is about 3.7 km beneath the volcano summit. The southern dyke reaches the depth of the modelled sill but is about 1 km south of the sill boundary.

## Complementary MT-INSAR results

To investigate Nyamulagira's inter-eruptive deformation behaviour with the goal to better constrain its magma plumbing system, we apply the 'StaMPS' technique to ENVISAT tracks that contain at least 20 images acquired between 2003 and 2010 over Nyamulagira.

Figure 20a shows the mean LOS velocity obtained for the ENVISAT, track 35, beam mode 2, descending orbit (incidence angle c. 18° from vertical) for the period January 2003 to March 2010. Most of the range increases (motion away from the satellite) correspond to cooling and compaction of historic lava flows. The most spectacular lava flow contraction is shown in Figure 20b, where an accumulation of Nyamulagira lava flows from eruptions in 1958, 1967, 1980 and 1991–1993 (Smets et al. 2010a) is still subsiding in a linear fashion with a LOS rate of about 25 mm a⁻¹.

The SE part of Nyamulagira's caldera (Fig. 20c) is similarly affected by persistent range increase behaviour, probably also owing to compaction of lavas. The shape of the deformation pattern is consistent with the shape of the 2001 lava flow inside the caldera; however, when looking closely at the time-series plot, one can identify c. 2 cm range

<table><caption>**Table 2.** Characteristics of interferograms used for modelling the 2010 Nyamulagira eruption, shown in Figure 17</caption><thead><tr><th>Satellite</th><th>Orbits</th><th>Time span (YYMMDD)</th><th>Incidence angle (deg)</th><th>H<sub>amb</sub> (m)</th><th>Number used in Figure 17</th></tr></thead><tbody><tr><td>ENVISAT</td><td>Ascending</td><td>091210–100218</td><td>31 (i5)</td><td>45</td><td>1</td></tr><tr><td>ENVISAT</td><td>Ascending</td><td>091229–100202</td><td>38 (i7)</td><td>430</td><td>2</td></tr><tr><td>ENVISAT</td><td>Descending</td><td>091210–100114</td><td>18 (i2)</td><td>43</td><td>3</td></tr><tr><td>ENVISAT</td><td>Descending</td><td>091207–100111</td><td>28 (i4)</td><td>74</td><td>4</td></tr><tr><td>ENVISAT</td><td>Descending</td><td>091223–100127</td><td>34 (i6)</td><td>3200</td><td>5</td></tr><tr><td>ENVISAT</td><td>Descending</td><td>091204–100108</td><td>38 (i7)</td><td>32</td><td>6</td></tr><tr><td>RADARSAT-2</td><td>Descending</td><td>091215–100201</td><td>35 (21)</td><td>97</td><td>8</td></tr><tr><td>ALOS</td><td>Ascending</td><td>091217–100201</td><td>34</td><td>704</td><td>9</td></tr></tbody></table>
<figure><img src="image_25.png" /><figcaption>**Fig. 19.** (a) Observed (left), modelled (middle) and residual (right) ENVISAT interferograms based on a model of two dyke intrusions and a deflating reservoir associated with the 2010 eruptive fissures (bold lines). Nyamulagira's summit caldera is outlined with a dashed white line. Note that all eight interferograms mentioned in the text and shown in Figure 17 were simultaneously inverted, but only two are shown here for clarity. (b) Geometry of the two dykes and reservoir model. Views are as in Figure 9.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Fig. 20.** Mean LOS velocity, MLV (in mm a⁻¹), from ENVISAT descending SAR data (incidence angle c. 18°) for the Nyamulagira area corrected for DEM and orbital errors and spanning 16 January 2003 to 25 March 2010. Lettered labels (b to f in part (a)) give the locations of time series shown in subsequent (b) to (f) figure parts. The reference area is indicated by the white rectangle. Red velocities correspond to range increase (ground moving away from satellite), and blue to range decrease. Nyamulagira's summit caldera is outlined in white. Red lines in time series plots indicate Nyamulagira eruptions and blue line marks the 2008 Bukavu earthquake.</figcaption></figure>

decreases just before the 2006 and 2010 eruptions. The small range decrease observed just before the 2010 eruption is confirmed by other ENVISAT time-series results (Wauthier 2011; Samsonov & d'Oreye 2012), and similar range decreases before

the 2006 and 2010 eruptions are also present NW of the caldera (Fig. 20e) and on the upper SE flank (Fig. 20d) along the area of weakness linking Nyiragongo and Nyamulagira volcanoes. The poor density of data early in the time series does not
allow us to interpret deformation behaviour before
the 2004 eruption.

A period of particularly high seismic activity
occurred during six months following a 5.9-
magnitude earthquake in Bukavu (128 km SSW of
Nyamulagira) on 3 February 2008 (d'Oreye et al.
2011). Small perturbations of around 2–3 cm in
the deformation behaviour can be observed in the
time-series results (Fig. 20) following the Bukavu
earthquake.

Discussion

Nyamulagira magma transport and budget

No magma was erupted in the caldera during flank eruptions in 1996 and 2006; however, the 2002, 2004 and 2010 eruptions were associated with erup- tive fissures and lava fountains on the caldera floor. Each of those three eruptions also included flank activity simultaneous with, or just after, the opening of the summit vents (on the north flank in 2002 and 2004 and on the SE flank in 2010). The flank eruptive fissure during the 2010 eruption was

located only c. 3 km to the north of the 2006
eruptive fissure.

Dyke overpressures found by modelling the 1996, 2002, 2006 and 2010 eruptions are between 0.2 and 3 MPa, with the northern dyke segment of the 2010 eruption defining the upper bound. Wadge & Burt (2011) suggested that dyke orien- tation might be controlled by tectonic stresses, with dykes orthogonal to the EAR requiring greater overpressures. Indeed, the fissure and dyke orientations for the 1996, 2002, 2006 and southern part of the 2010 eruptions are subparallel and could be favourably oriented with respect to the least principal stress or local rift extension (Wadge & Burt 2011). In such a model, the northern 2010 dyke might be less favourably oriented with respect to the least principal stress or local rift extension (Wadge & Burt 2011), and thus had a greater modelled overpressure.

The 1996 and 2002 eruptions, which both occurred on the north flank of the volcano, were not associated with any reservoir deflation. In contrast, the 2006 and 2010 eruptions, which both occurred on the SE flank of the volcano, were associated with a deformation of a shallow magma

<figure><img src="image_27.png" /><figcaption>**Fig. 21.** (a) NW–SE cross-section of topographic and best-fit source meshes inferred from the 1996, 2002 and 2010 modelling. (b) Conceptual sketch of Nyamulagira's magma plumbing system, adapted from (Pouclet 1976) and (Smets et al. 2010a). A–D show possible magma storage areas and transport pathways. Figure is not to scale.</figcaption></figure>
reservoir. Preliminary modelling of the 2006 InSAR data suggests that an inflating spherical or sill-like reservoir is present about 1.5 km below the summit, while in 2010 the summit reservoir is deflating and modelled at a depth of 3.7 km. Despite the difference in modelled depths, the results argue for magma storage at a few kilometres depth beneath Nyamulagira's summit (Fig. 21a). Seismic data (Hamaguchi & Zana 1983) suggest a magma reservoir between 3.5 and 7 km beneath Nyamulagira's summit. The surface projection of this inferred reservoir is consistent with the location of the depression in the SW part of Nyamulagira's caldera (Fig. 2). The magma budget between the reservoir and the erupted lava (Table 1) indicates that the reservoir is probably connected to a deeper source and is probably storing magma temporarily before eruptions.

According to our InSAR time series results, the caldera floor was affected by range decreases (motion towards the satellite) a few months before the 2006 and 2010 eruptions, while an overall range increase occurred during other periods (Fig. 20c, e). Thanks to a new multi-geometry time series method, Samsonov & d'Oreye (2012) identified unambiguous pre-eruptive deformation more than two weeks prior 2010 and also possibly for the 2004 and 2006 eruptions. Toombs & Wadge (2012) found 1–2 cm of inflation in a 10 km-wide zone centred on the caldera during the six months preceding the 2010 eruption. The pre-eruptive deformation could indicate addition of magma to a shallow reservoir a few months before each eruption. Mavonga et al. (2006, 2010) observed that LP and hybrid earthquakes preceded the 2004 and 2006 eruptions of Nyamulagira. About 10 months prior to each eruption, a steady increase in seismicity was observed, which is attributed to the movement of magma in a deep conduit. During the last one to two months before each eruption, the hypocentres of the LP earthquakes became shallower (Mavonga et al. 2006, 2010), which is consistent with supply to a shallow reservoir from a deep source. The largest range decrease was associated with the 2006 eruption. Overpressure in the magma reservoir was probably high, which could explain why the 2006 intrusion caused deformation over such a large area. The February 2008 Bukavu earthquake was also associated with a few centimetres of range decrease at Nyamulagira's summit area, suggesting that the earthquake perturbed the volcano's shallow plumbing system.

A schematic model for the plumbing system beneath Nyamulagira, consistent with the model proposed by Pouclet (1976), is shown in Figure 21b. In the model, area A is the region just beneath the volcano's summit, where intrusions feed summit eruptions and occasional lava lake activity

within the pit crater inside the caldera. Area B represents the main magma reservoir beneath the volcano. The deep reservoir is represented by area C, which feeds the shallower reservoir B. The depth of reservoir C is constrained by the seismic studies of Mavonga et al. (2010), who estimated the depth of this source to be 20–30 km. As indicated by graben subsidence on the SE flank of the volcano and by the time-series results indicating LOS subsidence of this region, magma could also be stored close to the surface in a dyke-like reservoir (area D in Fig. 21b). Such an elongated NW–SE magma storage zone is consistent with the NNW–SSE elongated reservoir inferred by Wadge & Burth (2011). Magma supply to a shallow reservoir (B and/or D) could be responsible for the range decrease observed before the 2006 and 2010 eruptions.

## Conclusions

A combination of 3D mixed boundary element method and inverse modelling was used to interpret InSAR displacements associated with eruptive activity in 1996, 2002, 2004, 2006 and 2010. All eruptions can be fitted by models incorporating dyke intrusions, and some (namely the 2006 and 2010 eruptions) require a magma reservoir beneath the summit caldera.

The northern dyke found in the modelling of the 2010 eruption has a larger overpressure than other Nyamulagira modelled dykes. This dyke might thus be less favourably oriented with respect to the local least principal stress and rift extension.

A magma reservoir beneath Nyamulagira has been identified from modelling deformation associated with the 2006 and 2010 eruptions. This reservoir is located 2–4 km beneath Nyamulagira's caldera, which is consistent with seismic data. It is probably connected to a deeper source and probably stores magma temporarily before eruptions.

The time-series results give us the following information: (1) the Bukavu 2008 earthquake could have induced a perturbation of the shallow plumbing system at Nyamulagira; and (2) the caldera and upper SE flank are subject to range increases in the months preceding the 2006 and 2010 eruptions, possibly indicating inflation of the shallow reservoir located beneath the caldera floor and SE flank.

We thus propose the following magma plumbing system at Nyamulagira: a deep reservoir (c. 25 km depth) feeds a shallower reservoir (c. 5 km depth). Proximal eruptions are fed by the shallow reservoir via dyke intrusions, while distal eruptions can be fed directly by the deep reservoir. We also infer the presence of a dyke-like conduit just beneath Nyamulagira's upper SE flank.
RADARSAT-1 data were provided by Canadian Space Agency via the Alaska Satellite Facility. RADARSAT-2 data covering Nyamulagira's 2010 eruption were provided by the Canadian Space Agency through SOAR project 5020. Other SAR data were made available through the ESA Cat-1 3224 and AO ALOS3690 projects, within the framework of the Belgian Science Policy (Belspo) and Luxembourg FNR, such as SAMAAV and GORISK (FNR/STEREOII/06/01). Precise orbits were provided by ESA. Interferograms were computed with the Repeat Orbit Interferometry PACkage (ROI-PAC) (JPL/ Caltech), DORIS (TU Delft) and GAMMA software. The SRTM digital elevation model was provided by USGS. Conventional interferograms were unwrapped using the SNAPHU algorithm. The authors thank M. Pritchard for his help concerning JERS and RADARSAT-1 data processing with ROI-PAC. We acknowledge Y. Fukushima and J.-L. Froger for the use of their inversion and meshing codes, respectively. Research by S.S. and B.S. was in part supported by the Luxembourg National Research Fund. We finally thank two anonymous reviewers and Editor J. Biggs for their thoughtful and useful comments.

## References

ALI, S. T. & FEIGL, K. L. 2012. A new strategy for estimating geophysical parameters from InSAR data: application to the Krafla central volcano in Iceland. *Geochemistry Geophysics Geosystems*, **13**, Q06005, http://dx.doi.org/10.1029/2012GC004112

AOKI, K. I. & YOSHIDA, T. 1983. Petrological and geochemical studies on the 1981–1982 lava from Nyamuragira volcano. In: *HAMAGUCHI, H. (ed.)* Volcanoes Nyiragongo and Nyamuragira: Geophysical Aspects. The Faculty of Science, Tohoku University, Sendai, 91–96.

AOKI, K., YOSHIDA, T., YUSA, K. & NAKAMURA, Y. 1985. Petrology and geochemistry of the Nyamuragira volcano, Zaire. *Journal of Volcanology and Geothermal Research*, **25**, 1–28.

BGVN 1996a. Nyamuragira, Bulletin Global Volcanism Network, Smithonian Institution, 21/10.

BGVN 1996b. Nyamuragira, Bulletin Global Volcanism Network, Smithonian Institution, 21/01.

BGVN 1998. Nyamuragira, Bulletin Global Volcanism Network, Smithonian Institution, 23/10.

BGVN 2000. Nyamuragira, Bulletin Global Volcanism Network, Smithonian Institution, 25/01.

BGVN 2001a. Nyamuragira, Bulletin Global Volcanism Network, Smithonian Institution, 26/01.

BGVN 2001b. Nyamuragira, Bulletin Global Volcanism Network, Smithonian Institution, 26/03.

BGVN 2001c. Nyamuragira, Bulletin Global Volcanism Network, Smithonian Institution, 26/12.

BGVN 2002a. Nyamuragira, Bulletin Global Volcanism Network, Smithonian Institution, 27/07.

BGVN 2002b. Nyamuragira, Bulletin Global Volcanism Network, Smithsonian Institution, 27/10.

BGVN 2003. Nyamuragira, Bulletin Global Volcanism Network, Smithonian Institution, 28/01.

BGVN 2004a. Nyamuragira, Bulletin Global Volcanism Network, Smithonian Institution, 29/04.

BGVN 2004b. Nyamuragira, Bulletin Global Volcanism Network, Smithonian Institution, 29/05.

BGVN 2006. Nyamuragira, Bulletin Global Volcanism Network, Smithonian Institution, 31/01.

BGVN 2007. Nyamuragira, Bulletin Global Volcanism Network, Smithsonian Institution, 32/01.

BÜRMANN, R., ROSEN, P. & FIELDING, E. 2000. Synthetic aperture radar interferometry to measure Earth's surface topography and its deformation. *Annual Reviews of Earth and Planetary Sciences*, **28**, 169–209.

BURT, M. L., WADGE, G. & SCOTT, W. A. 1994. Simple stochastic modelling of the eruption history of a basaltic volcano: Nyamuragira, Zaire. *Bulletin of Volcanology*, **56**, 87–97.

CALAIS, E., D'OREYE, N. *et al*. 2008. Strain accommodation by slow slip and dyking in a youthful continental rift, East Africa. *Nature*, **456**, 783–787, http://dx.doi.org/10.1038/nature07478

CAYOL, V. 1996. *Analyse élastostatique tridimensionnelle du champs de déformations des édifices volcaniques par éléments frontières mixtes*. PhD thesis, Université de Paris VII.

CAYOL, V. & CORNET, F. H. 1997. 3D mixed boundary elements for elastostatic deformations fields analysis. *International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts*, **34**, 275–287.

CAYOL, V., WAUTHIER, C., D'OREYE, N., KERVYN, F. & G. TEAM 2010. InSAR displacements associated with the November 2006 and January 2010 Nyamulagira eruptions. *AGU Fall Meeting*, Abstract #V32B-04.

CHAKRABATI, R., BASU, A. R., SANTO, A. P., TEDESCO, D. & VASELI, O. 2009. Isotopic and geochemical evidence for a heterogeneous mantle plume origin of the Virunga volcanics, Western rift, East African Rift system. *Chemical Geology*, **259**, 273–289.

CHENG, C. H. & JOHNSTON, D. H. 1981. Dynamic and static moduli. *Geophysical Research Letters*, **8**, 39–42.

COLCLOUGH, S. J. 2005. *Investigations of Nyamuragira and Nyiragongo volcanoes (DRC), using Interferometric Synthetic Aperture Radar*. PhD thesis, University of Cambridge.

DENAEYER, M. E. 1972. Les laves du fosse tectonique de l'Afrique Centrale (Kivu, Rwanda, Toro-Ankole). *Annales Musee Royal de Afrique Centrale Tervuren (Belgique), Series Geology*, **72**, 132.

D'OREYE, N., KERVYN, F. *et al*. 2008. Systematic InSAR monitoring of African active volcanic zones: what we have learned in two years, or an harvest beyond our expectations. *Proceedings of 2008 Second Workshop on USE of Remote Sensing Techniques for Monitoring Volcanoes and Seismogenic Areas*, Naples, 11–24 November 2008, http://dx.doi.org/10.1109/USEREST.2008.4740361

D'OREYE, N., GONZALEZ, P. J. *et al*. 2011. Source parameters of the 2008 Bukavu-Cyangugu earthquake estimated from InSAR and teleseismic data. *Geophysical Journal International*, **184**, 934–948, http://dx.doi.org/10.1111/j.1365-246X.2010.04899.x

EBINGER, C. J. 1989a. Tectonic development of the western branch of the East African Rift System. *Geological Society of America Bulletin*, **101**, 885–903.

EBINGER, C. J. 1989b. Geometric and Kinematic development of border faults and 539 accommodation zones, Kivu-Rusizi Rift, Africa. *Tectonics*, **8**, 117–133.

EBINGER, C. J. & FURMAN, T. 2002/2003. Geodynamical setting of the Virunga volcanic province, East Africa.
*In: Acta Vulcanologica, The January 2002 Eruption of Nyiragongo Volcano and the Socio-Economical Impact, 14 & 15.* Istituti editoriali e poligrafici internazionali, Pisa, 17–26.

FEIGL, K. L. & THURBER, C. H. 2009. A method for modeling radar interferograms without phase unwrapping: application to the M 5 Fawnskin, California earthquake of 1992 December 4. *Geophysical Journal International*, **176**, 491–504.

FOURNIER, T. J., PRITCHARD, M. E. & RIDDICK, S. N. 2010. Duration, magnitude, and frequency of subaerial volcano deformation events: new results from Latin America using InSAR and a global synthesis. *Geochemistry Geophysics Geosystems*, **11**, Q01003, http://dx.doi.org/10.1029/2009GC002558

FUKUSHIMA, Y. 2006. *Transferts de magma au volcan du Piton de la Fournaise déterminées par la modélisation 3D de données d'interférométrie radar entre 1998 et 2000.* PhD thesis, Laboratoire Magmas et Volcans, Université Blaise Pascal.

FUKUSHIMA, Y., CAYOL, V. & DURAND, P. 2005. Finding realistic dike models from interferometric synthetic aperture radar data: the February 2000 eruption at Piton de la Fournaise. *Journal of Geophysical Research*, **110**, B03206, http://dx.doi.org/10.1029/2004JB003268

FUKUSHIMA, Y., CAYOL, V., DURAND, P. & MASSONNET, D. 2010. Evolution of magma conduits during the 1998–2000 eruptions of Piton de la Fournaise volcano, Réunion Island. *Journal of Geophysical Research*, **115**, B10204, http://dx.doi.org/10.1029/2009JB007023

GRANDIN, R. A., SOCQUET, R. *et al.* 2009. September 2005 Manda Hararo-Dabbahu rifting event, Afar (Ethiopia): constraints provided by geodetic data. *Journal of Geophysical Research*, **114**, B08404, http://dx.doi.org/10.1029/2008JB005843

GUDMUNDSSON, A. & ANDREW, R. E. B. 2007. Mechanical interaction between active volcanoes in Iceland. *Geophysical Research Letters*, **34**, L10310, http://dx.doi.org/10.1029/2007GL029873

HAMAGUCHI, H. & ZANA, N. 1983. Introduction to Volcanoes Nyiragongo and Nyamuragira. *In: HAMAGUCHI, H.* (ed.) *Volcanoes Nyiragongo and Nyamuragira: Geophysical Aspects*. The Faculty of Science, Tohoku University, Sendai, 1–6.

HANSSEN, R. F. 2001. *Radar Interferometry: Data Interpretation and Error Analysis*. Kluwer, Dordrecht, Netherlands.

HEAD, E., SHAW, A. M., WALLACE, P. J., SIMS, W. W. & CARN, S. A. 2011. Insight into volatile behavior at Nyamuragira volcano (D.R. Congo, Africa) through olivine-hosted melt inclusions. *Geochemistry Geophysics Geosystems*, **12**, Q0AB11, http://dx.doi.org/10.1029/2011GC003699

HOOPER, A. 2008. A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches. *Geophysical Research Letters*, **35**, 302, http://dx.doi.org/10.1029/2008GL03465

HOOPER, A. 2010. Bayesian inversion of wrapped InSAR data for geophysical parameter estimation. *Proceedings of 'ESA Living Planet Symposium', Bergen, 28 June to 2 July 2010* (ESA SP-686, December 2010).

HOOPER, A., ÓFEIGSSON, B., SIGMUNDSSON, F., LUND, B., GEIRSSON, H., EINARSSON, P. & STURKELL, E. 2011. Increased capture of magma in the crust promoted by

ice-cap retreat in Iceland. *Nature Geoscience*, **4**, 783–786, http://dx.doi.org/10.1038/NGEO1269

HOOPER, A., SIGMUNDSSON, F. & PRATA, F. 2012. Remote sensing of volcanic hazards and their precursors. *Proceedings of the IEEE*, **100**, 2908–2930.

KAMPUNZU, M. M. A. B., SEBAGENZI, M. N., CARON, J. P. H. & VELLUTINI, P. 1982. Comparaison petrochimique entre les laves du champ Sud-Kivu (Bukavu) et du Nord-Kivu (Virunga), Zaire, *Annales de la Societe Geologique de Belgique*, **105**, 25–29.

KASAHARA, M. 1983. Near-field tilt measurements related to the 1981–1982 Nyamulagira eruption. *In: HAMAGUCHI, H.* (ed.) *Volcanoes Nyiragongo and Nyamuragira: Geophysical Aspects*. The Faculty of Science, Tohoku University, Sendai, 44–46.

KITAGAWA, S., FUKUI, K. & TAKAGI, A. 2007. The December 1996 eruption at Nyamuragira, Africa, viewed from JERS-1 InSAR. *In: Progress of Crustal and Deformation Research by Space Geodesy and Remote-sensing Techniques*, Kyoto University, 18–19 January 2007.

LU, Z., MASTERLARK, T. & DZURISIN, D. 2005. Interferometric synthetic aperture radar (InSAR) study of Okmok Volcano, Alaska, 1992–2003: magma supply dynamics and post-emplacement Lava flow deformation. *Journal of Geophysical Research*, **110**, B02403, http://dx.doi.org/10.1029/2004JB003148

MASSONNET, D. & FEIGL, K. 1995. Discrimination of geophysical phenomena in satellite radar interferograms. *Geophysical Research Letters*, **22**, 1537–1540.

MASSONNET, D. & FEIGL, K. L. 1998. Radar interferometry and its application to changes in the earth's surface. *Reviews of Geophysics*, **36**, 441–500.

MAVONGA, T. G. 2010. Crustal structure beneath two seismic broadband stations revealed from teleseismic P-wave receiver function analysis in the Virunga volcanic area, Western Rift. *Journal of African Earth Sciences*, **58**, 820–828.

MAVONGA, T., SADAKA, K. K., NYOMBO, L., OSODUNDU, E., JACQUES, D. N., ZANA, R. & DURRHEIM, R. J. 2006. Seismic activity prior to the may 8, 2004 eruption of volcano Nyamuragira, western rift valley of Africa. Studies of crustal structure, seismic precursors to volcanic eruptions and earthquake hazard in the eastern provinces of the Democratic Republic of Congo. *Journal of Volcanology and Geothermal Research*, **158**, 355–360.

MAVONGA, T., ZANA, N. & DURRHEIM, R. J. 2010. Studies of crustal structure, seismic precursors to volcanic eruptions and earthquake hazard in the eastern provinces of the Democratic Republic of Congo. *Journal of African Earth Sciences*, **58**, 623–633, http://dx.doi.org/10.1016/j.jafrearsci.2010.08.008

MISHINA, M., ZANA, N. & SAWA-SAWA, K. 1983. Geomagnetic Anomalies in the Virunga Volcanic Area and their Volcanological Implications. *In: HAMAGUCHI, H.* (ed.) *Volcanoes Nyiragongo and Nyamuragira: Geophysical Aspects*. The Faculty of Science, Tohoku University, Sendai, 47–54.

MOGI, K. 1958. Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. *Bulletin Earthquake Research Institute University of Tokyo*, **36**, 99–134.
OKADA, Y. 1985. Surface deformation due to a shear and tensile faults in a half-space. *Bulletin of the Seismological Society of America*, **75**, 1135-1154.

PARSONS, B. E., WRIGHT, T. J. *et al*. 2006. The 1994 Sefidabeh (eastern Iran) earthquakes revisited: new evidence from satellite radar interferometry and carbonate dating about the growth of an active fold above a blind thrust fault. *Geophysical Journal International*, **164**, 202–217.

POUCLET, A. 1976. Volcanologie du rift de l'Afrique centrale, le Nyamuragira dans les Virunga, essaie de magmatologie du rift. PhD thesis, University of Paris-Sud.

RUBIN, A. M. 1992. Dike-induced faulting and graben subsidence in volcanic rift zones. *Journal of Geophysical Research*, **97**, 1839–1858, http://dx.doi.org/10.1029/91JB02170

RUBIN, A. M. & POLLARD, D. D. 1988. Dike-induced faulting in rift zones of Iceland and Afar. *Geology*, **16**, 413–417.

SAMBRIDGE, M. 1999a. Geophysical inversion with a neighbourhood algorithm—I. Searching a parameter space. *Geophysical Journal International*, **138**, 479–494.

SAMBRIDGE, M. 1999b. Geophysical inversion with a neighbourhood algorithm—II. Appraising the ensemble. *Geophysical Journal International*, **138**, 727–746.

SAMSONOV, S. & D'OREYE, N. 2012. Multidimensional time series analysis of ground deformation from multiple InSAR data sets applied to Virunga Volcanic Province. *Geophysical Journal International*, **191**, 1095–1108.

SIMONS, M. & ROSEN, P. 2007. *Treatise on Geophysics, Interferometric Synthetic Aperture Radar Geodesy, Volume 3–Geodesy*. Elsevier, Amsterdam.

SMETS, B., WAUTHIER, C. & D'OREYE, N. 2010a. A new map of the lava flow field of Nyamulagira (D.R. Congo) from satellite imagery. *Journal of African Earth Sciences*, **58**, 778–786.

SMETS, B., KERVYN, M., KERVYN, F., D'OREYE, N. & WAUTHIER, C. 2010b. New insights into eruptive activity and lava flow hazard at Nyamulagira volcano, D.R.C., from a new GIS-based lava flow map. *AGU, 13–17 December 2010, San Francisco, CA, Abstract V11C-2288*.

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. *Geophysical Research Letters*, **35**, 6.

TARANTOLA, A. 1987. *Inverse Problem Theory: Methods for Data Fitting and Model Parameter Estimation*. Elsevier, Amsterdam.

TAZIEFF, H. 1977. An exceptional eruption: Mt. Nyiragongo, January 10th, 1977. *Bulletin of Volcanology*, **40**, 189–200, http://dx.doi.org/10.1007/BF02596999

TEDESCO, D., BADIALI, L. *et al*. 2007. Cooperation on Congo Volcanic and Environmental Risks, EOS, Transactions. *American Geophysical Union*, **88**, 177, 181.

TOOMBS, A. & WADGE, G. 2012. Co-eruptive and interruptive surface deformation measured by satellite radar interferometry at Nyamuragira volcano, D.R. Congo, 1996 to 2010. *Journal of Volcanology and Geothermal Research*, **245–246**, 98–122, http://dx.doi.org/10.1016/j.jvolgeores.2012.07.005

VERHOOGEN, J. 1948. *Les eruptions 1938–1940 du volcan Nyamuragira*, Institut des parcs nationaux du Congo Belge, Exploration du Parc National Albert, Missions J. Verhoogen (1938 et 1940), Fascicule 1.

WADGE, G. & BURT, L. 2011. Stress field control of eruption dynamics at a rift volcano: Nyamuragira, D.R.Congo. *Journal of Volcanology and Geothermal Research*, **207**, 1–15, http://dx.doi.org/10.1016/j.jvolgeores.2011.06.012

WAUTHIER, C. 2011. *InSAR Applied to the Study of Active Volcanic and Seismic Areas in Africa*. PhD thesis, University of Liege.

WAUTHIER, C., OYEN, A. M. *et al*. 2009. L-band and C-band studies of African volcanic areas. *Proceedings of 2009 Geoscience and Remote Sensing Symposium*, 2009 IEEE International, IGARSS, July 13–17 2009, Cape Town, 978-1-4244-3395-7/09, II 210–213, http://dx.doi.org/10.1109/IGARSS.2009.5418103

WAUTHIER, C., CAYOL, V., KERVYN, F. & D'OREYE, N. 2012. Magma sources involved in the 2002 Nyiragongo eruption, as inferred from an InSAR analysis. *Journal of Geophysical Research*, **117**, B05411, http://dx.doi.org/10.1029/2011JB008257

WRIGHT, T. J., EBINGER, C., BIGGS, J., AYELE, A., YIRGU, G., KEIR, D. & STORK, A. 2006. Magma-maintained rift segmentation at continental rupture in the 2005 Afar dyking episode. *Nature*, **442**, 291–294, http://dx.doi.org/10.1038/nature04978, PMID 16855588

YANG, X. M., DAVIS, P. M. & DIETERICH, J. H. 1988. Deformation from inflation of a dipping finite prolate spheroid in an elastic half-space as a model for volcanic stressing. *Journal of Geophysical Research* **93**, 4249–4257.

YUN, S., SEGALL, P. & ZEBKER, H. 2006. Constraints on magma chamber geometry at Sierra Negra Volcano, Galápagos Islands, based on InSAR observations. *Journal of Volcanology and Geothermal Research*, **150**, 232–243.