Nyamulagira’s magma plumbing system inferred from 15 years of InSAR

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 jeopardized 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 techniques, particularly interferometric synthetic aperture 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 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, 2006and 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/currentactivity/november-2011-eruption,2012) suggest a deformed area of more than 250 km 2 and an erupted volume of at least 81.5 × 10 6 m 3 .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 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 2 (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.

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; 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 N1108E, with an estimated extension rate of 2.8 mm a 21 (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 are labelled.The white rectangle gives the extent of Figure 20.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 subcentimetre accuracy (Massonet & 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 (l 5.6 cm) interferograms having small temporal and perpendicular baselines (Fig. 4;d'Oreye et al. 2008), or L-band (l 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 2 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 2 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  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(Fukushima et al. , 2010)).Three parameters -Dip (dip angle), Shear (angle between the dip direction and the line connecting the  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 e (measured from horizontal; a horizontal ellipse will thus have a Dip e ¼ 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.Tractions are zero at the ground surface and are equal to overpressures in  dykes or reservoirs (Cayol & Cornet 1997).With the goal to constrain the stress field in the lower ductile 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 23 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 Nyamulagira, 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, constant dyke overpressures are therefore assumed.This boundary condition leads to a non-uniform distribution of displacements.The source geometries are inverted simultaneously with the tractions, leading 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 horizontal dimension of the deformation source (Cayol 1996).The density of the topographic mesh is fine close to the deformation sources, where displacement 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 eruption, 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 composed 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 deformation, leading to an underestimation of source overpressure (Fukushima et al. 2005).A good compromise between calculation time and model accuracy 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 simultaneously invert for the dyke geometry, location and uniform overpressure.As surface displacements are a non-linear function of the source parameters, 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 discrepancy between observed and 3D-MBEM modelled displacements (Fukushima et al. 2005): where C d is the data covariance matrix and u o and u c , formed by concatenating data from the different viewing geometries, are vectors of observed and modelled line-of-sight (LOS) displacements, respectively.
To make the misfit computation numerically manageable, u o is constructed by subsampling the unwrapped LOS displacements (Fig. 7b).NYAMULAGIRA'S MAGMA PLUMBING SYSTEM 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: where r is the distance between two pixels, s (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 24 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: 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 6 m 3 (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 preeruptive (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 2 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 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 × 10 6 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 2002InSAR in , 2004InSAR in , 2006   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 + 15.8 × 10 6 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 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(BGVN , 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 6 m 3 (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 longperiod (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 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 bestfitting model is a dyke that is longer at depth than the eruptive fissure, dipping east at shallow levels, and dipping 408 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 × 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 + 22.8 × 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.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(BGVN , 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 118 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 + 14.7 × 10 6 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 2 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 2 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 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 + 17.9 × 10 6 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   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 808 eastward, has a volume of 1.1 × 10 6 m 3 , and reaches a depth of 2.6 km beneath the summit of the volcano, while the southern dyke dips 708 westward, has a volume of 2.7 × 10 6 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).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.tions in 1958tions in , 1967tions in , 1980tions in and 1991tions in -1993tions in (Smets et al. 2010a) is still subsiding in a linear fashion with a LOS rate of about 25 mm a 21 .
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    allow us to interpret deformation behaviour before the 2004 eruption.
A period of particularly high seismic activity occurred during six months following a 5.9magnitude 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.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 orientation 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.

Nyamulagira magma transport and budget
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 NYAMULAGIRA'S MAGMA PLUMBING SYSTEM 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. (2006Mavonga et al. ( , 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(Mavonga et al. , 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.

Fig. 3 .
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.

Fig. 2 .
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).

Fig. 4 .
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.

Fig. 5 .
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).

Fig. 6 .
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).

Fig. 7 .
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.
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 eruptionsThe 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 + 23.1 × 10 6 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

Fig. 8 .
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.

Fig. 9 .
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.

Fig. 12 .
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.

Fig. 13 .
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. 188, H amb ¼ 49 m).(b) Ascending RADARSAT-1 interferogram (26 March 2004 to 30 June 2004, incidence angle c. 448, H amb ¼ 26 m).(c) Descending RADARSAT-1 interferogram (16 March 2004 to 7 August 2004, incidence angle c. 358, H amb ¼ 30 m).Nyamulagira's summit caldera is outlined with a dashed white line.

Fig. 14 .
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.

Fig. 15 .
Fig. 15.ENVISAT ascending interferogram (31 October 2006 to 5 December 2006, incidence angle c. 438, 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.

Fig. 16 .
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.

Fig. 17 .
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.SeeTable 1 for data description.

Figure
Figure 20a shows the mean LOS velocity obtained for the ENVISAT, track 35, beam mode 2, descending orbit (incidence angle c. 188 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 erup-tions in 1958tions in  , 1967tions in  , 1980tions in   and 1991tions in  -1993tions in   (Smets  et al. 2010a) is still subsiding in a linear fashion with a LOS rate of about 25 mm a 21 .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

Fig. 18 .
Fig. 18.(a) ENVISAT descending wrapped interferogram spanning 10 December 2009 to 14 January 2010, i7 beam mode (incidence angle c. 388), and (b) ascending wrapped interferogram spanning 29 December 2009 to 2 February 2010, i7 beam mode (incidence angle c. 378); 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.

Fig. 19 .
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.

Fig. 20 .
Fig. 20.Mean LOS velocity, MLV (in mm a 21 ), from ENVISAT descending SAR data (incidence angle c. 188) 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.
No magma was erupted in the caldera during flank eruptions in 1996 and 2006; however, the 2002, 2004 and 2010 eruptions were associated with eruptive 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

Fig. 21 .
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.
Table 1 for data description.

Table 2 .
Characteristics of interferograms used for modelling the 2010 Nyamulagira eruption, shown in Figure17