Magma sources involved in the 2002 Nyiragongo eruption, as inferred from an InSAR analysis

2.1. Geologic Setting

[6] Nyiragongo volcano (Figure 1) lies on the East African Rift (EAR). The 3000 km-long EAR has propagated from north to south, as shown by systematic variations in geophysical, geochemical and geomorphological data [Chorowicz, 2005]. Spreading rates decrease from 6.5 mm/year in the north of the rift to 1.5 mm/year in the south of the rift [Stamps et al., 2008]. In the Virunga Volcanic Province (VVP), where the rift trend is NNE, with localized NE trends, the spreading rate has an intermediate value of 2.8 cm. The rift in the VVP is most probably a half-graben with a marked normal fault on the western side and no clear normal bounding fault on the eastern side (Figure 1) [Ebinger, 1989a].

[7] Nyiragongo and Nyamulagira volcanoes, ∼15 km apart, are the youngest and most westerly volcanoes of the VVP (Figure 1). Distribution of cones and eruptive fissures at Nyamulagira and Nyiragongo generally follows the NE rift trend, but there is also a NNW-trending fracture network (Figure 1), which may be due to the stress field generated by the interplay between the magmatic systems of both volcanoes [Gudmundsson and Andrew, 2007] or to the reactivation of Precambrian basement faults [Smets et al., 2010b]. Nyamulagira and Nyiragongo volcanoes are located at the intersection between these two trends. The location of the volcanic centers with respect to the Kivu border faults indicates interaction between the volcanism and faulting [Ebinger, 1989b].

[8] While Nyamulagira erupts every 2–4 years [Smets et al., 2010a], only two historical eruptions of Nyiragongo have been recorded; in 1977 [Tazieff, 1977] and in 2002 [Komorowski et al., 2002/2003; Tedesco et al., 2007]. Nyiragongo is one of the few volcanoes on Earth to host an active lava lake. In both eruptions, the summit lava lake drained and fracture networks opened on the southern flank of the volcano (Figure 2).

2.2. Magma Reservoirs at Nyiragongo

[9] Magma reservoirs beneath Nyiragongo are still poorly constrained. Demant et al. [1994] explain the petrography of the lavas based on two reservoirs, located at different depths, and periodically connected. The shallow reservoir, where the final fractional crystallization occurs, is assumed to be permanently connected to the active lava lake. Estimates of the shallow reservoir depth, from leucite crystallization pressures, place it at a depth varying from 4 km [Platz et al., 2004] to less than 1 km [Louaradi, 1994]. The depth of the deep reservoir, where the initial crystallization occurs, is estimated at 10–14 km, based on petrographical analysis of fluid and melt inclusions in olivine and clinopyroxene phenocrysts [Demant et al., 1994]. This location is consistent with seismic studies which attribute the a-seismic region, which lies between the surface and a depth of 10–14 km, to a magma complex [Tanaka, 1983]. The existence of at least two magma reservoirs is confirmed by a seismic study by Shuler and Ekström [2009], which showed that a series of five unusual earthquakes of moderate size recorded during the 2002 eruption could be attributed to the successive collapse of the roofs of the two connected magma reservoirs.

2.3. The January 2002 Nyiragongo Eruption

[10] On 17 January 2002, Nyiragongo erupted along a 20 km-long fracture network extending southward from the northwestern flank of the volcano, toward the airport and city of Goma. Part of the 1977 fracture network was reactivated (Figure 2). A small radial fissure opened on the northwestern flank of the volcano, while a 12 km-long N-S fracture system opened on the southern flank, almost parallel to the main rift-related faults (Figure 2). The observed opening is about 1–3 m. The eruption lasted for 12 [Tedesco et al., 2007] to 48 h [Komorowski et al., 2002/2003]. The erupted volume is estimated to be in the order of 14–34 × 10⁶ m³ [Tedesco et al., 2007].

2.4. Available in Situ Data

[11] Field surveys, conducted after the eruption, provided lava samples from different parts of the eruptive fracture system. Based on ²³⁸U and ²³²Th series radioactive disequilibria data, Tedesco et al. [2007] inferred that up to three different magma sources were involved in the eruption. The northern vents emitted degassed lava drained from the lava lake, possibly combined with input from a shallow magmatic reservoir located below the main Nyiragongo edifice. The southern part of the eruptive fissure, which opened at the end of the eruption, released gas-rich lavas as indicated by the lava fountains. This observation, together with the disequilibria analysis, indicates that this magma originated from a deeper source, probably located beneath the city of Goma.

[12] After the eruption, carbon dioxide and methane emanations were observed around the eruptive fissures and in several locations in Goma between 300 and 800 m from the lava flows [Komorowski et al., 2002/2003]. Methane bursts were also reported in Goma, south of Munigi (Figure 2) during the eruption [Komorowski et al., 2002/2003]. This methane could be stored in pockets [Tedesco et al., 2007] located within the sediments filling the Lake Kivu basin and might be released continuously via tectonic fractures throughout the whole region [Komorowski et al., 2002/2003]. The sudden fracturing associated with the 2002 eruption could have triggered the degassing.

[13] The two operating seismic stations showed that precursory seismicity started ten months before the eruption [Kavotha et al., 2002/2003]. It was characterized by the onset of continuous volcanic tremor, and an increased number of long-period and short-period earthquakes, sometimes followed by increased fumarolic activity [Tedesco et al., 2007]. Magma might thus have been supplied to a reservoir in the months preceding the eruption [Kavotha et al., 2002/2003]. Between 4 and 16 January 2002, and after the effusive activity had stopped, seismic activity remained high [Kavotha et al., 2002/2003]. About 100 tectonic earthquakes with magnitudes greater than 3.5 were recorded in the five days following the eruption over the whole Lake Kivu area [Tedesco et al., 2007], culminating in a 6.2 M_w event on 24 October, 2002 [Mavonga, 2007]. Although the seismic network was not capable of accurately determining the location of the seismic events, they generally cluster along the faults on the western side of the Kivu basin, as found for previous recorded events [Zana and Hamaguchi, 1978; Mavonga, 2007]. The Harvard Centroid Moment Tensor (CMT) focal mechanisms for the three largest January 2002 earthquakes (>5 M_w) indicate almost pure normal faulting mechanisms.

[14] Building damage was reported during the two weeks after the eruption, especially in Gisenyi (Figure 2b) [Komorowski et al., 2002/2003]. During field investigations carried out in January 2010, the inhabitants of the damaged houses indicated that the initiation of the cracks was contemporaneous with the January 2002 Nyiragongo eruption. New cracks also opened up several years after the eruption, some appearing in early 2009. The Rwandese Red Cross, which has been measuring the extension of certain cracks in the Gisenyi Hospital, confirmed that crack widening is still ongoing. In some cases, vertical movement was noted, with the western side moving down.

[15] After the eruption, local inhabitants reported a noticeable rise in the water level of the lake along the shoreline in the Goma area. Eleven days after the eruption, water level changes were estimated by measuring the depth of the immerged algae zone (∼30 cm thick) that develops just below the water surface (J. Durieux, personal communication, 2005). These measurements were repeated at irregular intervals until March 2002 [Komorowski et al., 2002/2003; Tedesco et al., 2007] and showed an asymmetric increase in the water level of Lake Kivu, which reached a maximum of about 37 cm at the Goma shore (Profile A in Figure 3). Due to the agitation of the water surface by waves, and given possible additional sources of perturbations such as water discharge at the 100 km-distant Ruzizi hydroelectric dam, seasonal effects, seiches, etc., uncertainties can conservatively be estimated to be of the order of 30 cm. This water level change has been attributed [Komorowski et al., 2002/2003] to subsidence of the rift system around Goma as the asymmetry is consistent with the asymmetry of the graben in which the VVP is located [Ebinger, 1989a].

4.1. Mechanical Properties of the Crust

[28] The edifice is assumed to be linearly elastic, homogeneous and isotropic. However, elastic properties vary over the rift depth and width. Although there is no data to document the lateral variation of the elastic properties, a seismic velocity profile of the Virunga area [Mavonga, 2010] showed that seismic velocities increase from 3.5 km/s at sea level to 7.3 km/s at 7.5 km b.s.l. (Figure 6a). Experimental values of Poisson's ratio, ν, for basalts range from 0.22 to 0.28 for low to high confining pressure, corresponding to a depth of several km [Birch, 1966]. We take an average Poisson's ratio of 0.25. The Young's modulus, density and Poisson's ratio contribute to the seismic velocities. Young's modulus can be determined from:

where V_p is P wave the seismic velocity, E_d is the dynamically determined Young's modulus and ρ_r is the rock density. Based on studies in Hawaii [Zucca et al., 1982], we assume that the density varies from 2500 kg/m³ at shallow depth, where vesiculated lava flows are present, to 2900 kg/m³ at greater depth where the rocks are compacted (Figure 6b). Density at shallow depth is consistent with the value of 2500 kg/m³ estimated from a gravity survey around Nyiragongo and Nyamulagira volcanoes [Mishina, 1983]. In the absence of data we assume the density variation to be linear with depth.

[29] From seismic velocities, the dynamic Young's modulus is estimated to increase with depth from 23 GPa to 124 GPa (Figure 6c). Laboratory experiments [Cheng and Johnston, 1981; Adelinet et al., 2010] show that statically determined values of Young's modulus, E_s, are always less than dynamically determined values, E_d, and that the E_s/E_d ratio increases with confining pressure, and decreases with porosity and crack density. Here, confining pressure is estimated by integrating the product of the density and the depth (Figure 6d).

[30] We further consider an elevated porosity and crack density close to the surface, leading to a lower ratio of E_s/E_d of 0.1 which was determined by Cheng and Johnston [1981] for tuffs at very low confining pressures. The E_s/E_d ratio for larger confining pressures is estimated from the study of Adelinet et al. [2010]. It is assumed to increase from E_s/E_d = 0.25 at sea level to 0.8 at a confining pressure of 200 MPa (Figure 6e). Thus, our first order estimate is that the Young's modulus varies from 2.3 GPa close to the surface to 100 GPa at 9.5 km below the surface (Figure 6f).

[31] Preliminary models [Wauthier, 2007] showed that the dike associated with the 2002 eruptive fissure is shallow. As our modeling method is for a homogeneous medium, we must consider a single Young's modulus value in the computations. We take a value of 5 GPa, corresponding to the Young's modulus value in the upper two kilometers of the crust, which is the average depth of the presumed eruption dike.

4.2. Three-Dimensional Boundary Element Modeling

[32] The numerical modeling method is a mixed boundary element method (MBEM) [Cayol and Cornet, 1997]. Boundary element methods are used to solve linear partial differential equations that have been formulated as integral equations. Linear integral equations relating Green's function to prescribed boundary conditions are obtained and simultaneously solved for the unknown boundary conditions, so that mechanical interactions between sources are taken into account. The MBEM computes the three-dimensional displacements caused by sources such as tensile cracks, faults or reservoirs of arbitrary geometry located beneath a realistic topography. Nyiragongo volcano has steep slopes of more than 20 degrees. At Piton de la Fournaise volcano, where slopes are similar to Nyiragongo, Fukushima et al. [2005] found that neglecting the topography resulted in an overestimation of the overpressure (30%), height (30%) and volume (20%) of the modeled dike intrusion.

[33] Boundaries are meshed by planar triangular elements. The topographic mesh is constructed from the DEM used for the InSAR processing and also includes a rough bathymetry [Capart, 1955; Schmid et al., 2010]. The size of the mesh is about five times greater (200 km in radius) than the deformation extent, so that the limited extension of the ground surface has only a small influence (around 5%) on the computed displacements [Cayol, 1996]. The topographic mesh is dense close to the eruptive fissures in the Goma area where displacement gradients are greatest, and it becomes coarser further away (Figure S2 in the auxiliary material).

[34] Cayol and Cornet [1997] and Fukushima et al. [2005] showed that coarse meshing of the sources generally induces an overestimation of the amplitude of surface deformation, leading to an underestimation of source overpressure. We determine mesh intervals for the source by progressively diminishing the source mesh interval so that it no longer influences the results. A 400 m mesh interval is used for the dike associated with the 2002 eruptive fissure and elliptic reservoirs, while an 800 m mesh interval is used for the buried dike and faults.

[35] When a dike is connected to the ground surface, it's shape is assumed to be close to that of a quadrangle. The top part of the dike corresponds to a single smooth curve through the segmented eruptive fissure. The bottom part of the dike corresponds to a straight line defined by six parameters Dip, Shear, Botelv, Botlen, Twist and Botang (Figure 7a). These parameters are chosen in order to be able to restrict the model search to mechanically plausible models [Fukushima et al., 2005]. We do not consider the individual fissure segments mapped after the eruption (Figure 2b). Indeed, according to St Venant's principle [Love, 1927], they only influence displacements close to the eruptive fissure, where InSAR data are incoherent.

[36] To keep the number of parameters reasonably low for an efficient inversion, buried fractures are assumed to have a simpler quadrangular shape, with straight top and bottom sides (Figure 7b), leading to seven geometrical parameters, similar to the frequently used rectangular dislocation model of Okada [1985], which has three parameters when the top of the quadrangle is fixed and seven parameters when it is unknown.

[37] The prescribed boundary conditions are stress vectors that are assumed to correspond to perturbations of an initial state of stress. Indeed, a previous study [Zeller and Pollard, 1992] showed that boundary element solutions obtained with stress boundary conditions were physically more realistic than those obtained with displacement boundary conditions. It is assumed that there is no stress change at the ground surface, so that stress vectors are zero at the surface. Dikes are assumed to open as a response to magma overpressure, which is defined as being the difference between the magma pressure and the normal stress exerted by host rocks on the dike surface. Faults are assumed to slip as a response to drops in shear stress [Cayol and Cornet, 1997].

[38] We assume that each dike is submitted to a uniform overpressure. Since dikes are fluid-filled fractures, all parts of a dike are hydraulically connected. We further assume that the host rock stress is homogeneous close to the dike, so that the overpressure is uniform or varies linearly with depth. A depth-dependent gradient could have been assumed, but a previous study [Fukushima et al., 2005] showed that the range of acceptable gradients was so great that it could not be constrained. Moreover, we do not have a good a priori knowledge of the host rock stress, thus, of the overpressure. Using a constant stress model for dikes leads to a smoothly varying opening on the fracture [Cayol and Cornet, 1997], consistent with results obtained from linear kinematic inversions of opening distributions [Yun et al., 2006; Pallister et al., 2010]. We invert the dike and fault geometries simultaneously with uniform stress vectors, leading to a minimum of 7 parameters per structure.

4.3. Nonlinear Inversions

[39] Following Fukushima et al. [2005], we use a near neighborhood inversion algorithm that consists of two stages: search and appraisal. In the search stage [Sambridge, 1999a], the model space is explored by iteratively sampling regions where data are well fitted. Multiple minima are determined. A misfit function is used which quantifies the discrepancy between observed and modeled displacements. It is written as:

where u_o and u_m are vectors of observed and modeled line-of-sight (LOS) displacements, respectively, formed by concatenating data from the different beams, and C_d is the data covariance matrix.

[40] To make the misfit computation numerically manageable, vectors of observed and modeled displacements u_o and u_m are constructed by subsampling the observed and modeled LOS displacements at the same points. Data are subsampled at circular gridded points in such a way that the subsampling is dense (the average distance between points is 800 m) in areas where displacement gradients are steep, i.e., close to the 2002 eruptive fissure and in the vicinity of Goma, and coarser further away. The choice of a subsampling circular grid method was guided by the results of Fukushima et al. [2005], who showed that, for dikes, the subsampling method does not significantly affect the inversion result. Each pixel with a coherence threshold of less than 0.15 is rejected from the subsampling process (Figure S3 in the auxiliary material).

[41] The vector of observed LOS displacements is such that u_oⁱ = u_oⁱ′ + offⁱ, where u_oⁱ′ and u_oⁱ are the observed displacements of the ith interferogram before and after it is shifted by the amount offⁱ. The shift of each of the interferograms offⁱ is determined for each forward model calculated during a near-neighbor inversion, by minimizing the misfit (equation (2)) between observed and modeled displacements ∂χ²/∂offⁱ = 0. As it is a linear model parameter, the shift is given by:

where is the inverse of the Covariance matrix of displacements in the ith interferogram, and eⁱ = (1..1)^T is a unity vector with the same dimension as u_oⁱ′ and u_mⁱ.

[42] InSAR data contain spatially correlated phase noise caused by the random fluctuation of atmospheric phase delays from one SAR image to another. Fukushima et al. [2005] found they can be appropriately modeled by an exponential autocorrelation or covariance function, expressed as

where C is either the autocorrelation or covariance function, r is the distance between two pixels, σ_d² is the variance of the noise, and a is the correlation distance [Fukushima et al., 2005]. Following Fukushima et al. [2005], we compute the autocorrelation function for the undeformed area of each interferogram and estimate the variance and correlation distance corresponding to the best exponential curve fit for each interferogram (Table 1). We then use these values to compute the full covariance matrix C_d which we use for the misfit computation (2).

[43] The appraisal stage [Sambridge, 1999b] involves calculations of posterior probability density functions using misfit values calculated during the search stage. It allows for the estimation of several model statistical characteristics such as the mean model and marginal posterior probability density functions (PPDs). The one-dimensional marginal PPDs provide confidence intervals for the model parameters, while the two-dimensional marginal PPDs indicate trade-offs between pairs of model parameters. The applicability of the method has been confirmed by synthetic tests [Fukushima et al., 2005].

5.1. Rejected Models

[47] Models considering one single dike (Table 2, column 1, and descriptions, and Figure S4 and Table S1 in the auxiliary material) or a dike and a deflating reservoir beneath the Goma area (Table 2, column 2, and descriptions, and Figure S5 and Table S2 in the auxiliary material) do not fit the InSAR data satisfactorily, leading to RMS errors of 4.8 and 4.2 cm, respectively. Thus they were rejected.

[48] Models considering a dike and several faults were also rejected (Figure S6 and Table S3 in the auxiliary material). In contrast to the 2007 Natron basin seismo-magmatic event [Baer et al., 2008; Calais et al., 2008; Biggs et al., 2009], field observations performed after the Nyiragongo eruption [Komorowski et al., 2002/2003] showed no significant fault trace, with the exception of the two 2.5 km-long 1/3 fringe discontinuities visible in the RADARSAT-1 descending interferogram in the Gisenyi area (Figure 5 and 17). Inversions considering a dike associated with the 2002 eruptive fissure and one to three westward-dipping normal faults located to the west of this dike led to a poor data fit with a RMS error of 3.5 cm (Table 2, column 3).

5.2. Preferred Model

[49] Our best fit model consists of a shallow dike associated with the eruptive fissure, and a second, deeper dike, with its northern end located beneath the Nyiragongo edifice, extending below Lake Kivu, and ending about 15 km from Idjiwi Island. Because of the lack of InSAR measurements over the lake surface, we assessed the position of the southern end of the deep dike based on the coherent InSAR data along the shore. It was fixed by trials and errors in the preliminary inversions. A reasonable data fit was found with the southern end of the deep dike located at 743 N and 9,805 E km UTM (Table 3b). The northern end of the deep dike is defined by two parameters (Strike and Toplen) that are inverted.

[50] In order to determine the most likely model, we performed six inversions, increasing the number of inverted parameters from 9 to 14 (Table 2, models (a) to (f)). Each inversion involved several thousand models (Table 2). The more complex the model, the larger the number of forward models required for the inversion to converge (Table 2). The 14 parameter model (Figure 8 and Table 2, model (f)*) is the most likely as it has the lowest AIC. It has the best data fit as it is associated with the lowest RMS error (2.3 cm, Table 2). With this model, the dike connected to the eruptive fissure is subvertical and shallow, with an average height of 2 km and a sub-horizontal bottom, making the dike nearer to the surface beneath Nyiragongo volcano (Figure 8 and Table 3a) than beneath the city of Goma. The dike is associated with an overpressure of 0.9 MPa. The opening distribution corresponding to this overpressure is maximum at the ground surface (Figure 8b), where the mean opening is 3 m, a value which is consistent with field measurements [Komorowski et al., 2002/2003]. The other dike is deeper and has a 20 km-long top at a depth of 3.3 km below the surface, an average height of 6 km, and a 40 km-long bottom (Figure 8a and Table 3b). Like the shallower dike, the deeper dike has a low overpressure of 0.4 MPa. The corresponding opening is at its maximum toward the center of the dike, reaching 1 m, with an average value of 0.7 m (Figure 8b). The two-dimensional PPDs (Figure S7 in the auxiliary material) show that there are no significant trade-offs between the model parameters.

[51] To the west of the 2002 eruptive fissure, residuals (Figure 9) are lowest on the shortest time-span interferogram. This indicates that time-dependent processes, such as a slow motion of the rift border faults or magma movements associated with Nyamulagira's 2002 eruption, might also have taken place. Other residuals, in the Goma area and east of Lake Kivu, are probably due to the oversimplified shapes of the dikes considered in the inversions.

6.1. Bias Induced by the Homogeneous Elasticity Hypothesis

[52] A low value of the Young's modulus was considered for the inversions of both dikes (E = 5 GPa), corresponding to the Young's Modulus in the upper two kilometers of the crust. Although this is a reasonable assumption for the shallow dike, such a low value could lead to an underestimation of the depth of the deeper dike. Models with depth-dependent elastic properties, with more compliant layers at shallow depth and stiffer layers at greater depth, always lead to over-estimations of source depths, whether these sources are shear faults [Arnadottir et al., 1991; Montgomery-Brown et al., 2009] or magmatic reservoirs [Masterlark, 2007].

[53] This low Young's modulus might also lead to an underestimation of the overpressure of the deep dike. Considering the depth-dependent rigidity estimated earlier for this area (8 and Figure 6f), the Young's modulus reaches an average value of 90 GPa at 4 km depth, which is the mean depth of the deep dike (Table 3b). When inverting displacements for overpressures in linear elastic media, the determined overpressures are linear functions of the Young's modulus. Thus, if the Young's modulus is E = 90 GPa instead of E = 5 GPa used in the inversions, the overpressure of the deep dike is 7.4 MPa, instead of 0.4 MPa.

6.2. Consistency With Lake Level Data

[54] In order to compare lake level data with the InSAR observations and our best fit deformation model, the InSAR displacements from the three beams were simultaneously used to compute vertical displacements along the lake shore following the procedure of Wright et al. [2004]. Before computing vertical displacements the InSAR displacements were shifted by the amount that minimized the misfit between the observed data and the best fit model (equation (3)). Standard deviations and unit LOS vectors considered for the displacement calculations are given in Table 1. Comparison of the vertical displacement of the lake shore as determined from the best fit model, the InSAR data and the water level measurements (Figure 3) shows a similar overall shape of subsidence characterized by a maximum subsidence at approximately the same location. Lake level data show a maximum subsidence of 32 cm relative to the easternmost point of the northern shore of Lake Kivu [Tedesco et al., 2007], greater than the modeled (8 cm) or InSAR measured (5 cm) relative subsidences. In addition, whereas the lake level data indicate only broad subsidence, the InSAR data and associated model predict uplift on each side of the subsiding sector. Due to the large uncertainties in the lake level measurements, we attribute the discrepancies between the model, InSAR data and the lake level measurements to errors in the latter.

6.3. Discontinuities in the Gisenyi Area

[55] In order to test whether these faults could have been activated by the 2002 dikes, we have computed the expected static stress changes as follows. Coulomb stress changes ΔS are given by

where Δτ is the shear stress change on failure planes with a given orientation and rake (positive in the sense of slip), Δσ is the change in normal stress, and μ′ is the effective friction coefficient which incorporates the effect of pore fluids. Following Stein et al. [1992] and King et al. [1994] a value of 0.4 is assumed. No attempt is made to investigate the influence of this parameter as several studies have concluded that Coulomb stress modeling is only marginally dependent on the effective friction coefficient values [King et al., 1994; Parsons et al., 2012]. As the Gisenyi discontinuities are parallel to the normal faults at the eastern rift border (Figure 2), which are west dipping, they are assumed to be the surface expression of the border faults. A dip of 60° is assumed as it is the most common dip for normal faults. A Coulomb stress calculation performed with the MBEM (Figure 10) shows that the discontinuities are located in an area where slip on normal faults would be promoted. These faults are too far from the dikes to be located in the area of greatest Coulomb stress increase. Thus, most probably, these faults were not created by the dike intrusions, but reactivated by them. As the observed displacements have a small amplitude and extend over a small area, it is likely that the faults were only reactivated to a shallow depth.

6.4. Magma Origin

[56] The two-dike model with a shallow and deep dike is consistent with the two reservoirs inferred from petrographical and seismic analysis. The shallow dike goes down to 3 km below the summit, a depth consistent with the shallow reservoir depth estimated to lie between 4 km [Platz et al., 2004] and 1 km [Louaradi, 1994] based on leucite crystallization pressure. The 12 km maximum depth below the surface for the deep dike is consistent with the calculated value for the deep reservoir of 10–14 km depth, as determined from petrographical analysis on fluid and melt inclusions in olivine and clinopyroxene phenocrysts [Demant et al., 1994], as well as the a-seismic region [Tanaka, 1983].

[57] The model is also consistent with ³⁸U and ²³²Th series radioactive disequilibria data, calculated from lavas emitted during the eruption [Tedesco et al., 2007], which suggest that the eruption was fed by up to three different magmatic sources. These authors infer that the northern lava flows were probably supplied by the lava lake and a shallow magma reservoir, while the southern lava flows were related to a second, deeper reservoir, probably located beneath Goma city and Lake Kivu. The shallow dike volume is consistent with this model as the inverted shallow dike has a slightly greater volume (74 × 10⁶ m³, Table 4) than the volume of lava that drained from the summit crater (60 × 10⁶ m³, Table 4 [Komorowski et al., 2002/2003]). Thus, it could have been supplied by other sources, such as a shallow or a deep reservoir.

[58] We do not see any evidence of a shallow or deep reservoir deflation, as suggested by the unusual earthquakes recorded during the eruption [Shuler and Ekström, 2009], but these deflations could have been hidden by the large displacements associated with the dike intrusion, or by the lack of coherent InSAR data on the steep vegetated slopes of the volcano and on Lake Kivu. The pressure in the reservoir could also have been restored by gas exsolution in the magma chamber [Rivalta and Segall, 2008].

[59] In our preferred model, the deep dike is not connected to the ground or to the other dike (Figure 8). However, the two dikes are only separated by about 2 km vertically, so the presence of a magma conduit between them is highly likely. The connection between both dikes is probably via another dike, as this is the most efficient means of transporting magma through the lithosphere [Rubin, 1995; Delaney and Pollard, 1981]. This dike might be too narrow to cause detectable InSAR displacements. Indeed, the data fit breaks down when the southernmost part of the shallow dike is shortened by 5 km and replaced with a 5 km-long dike connecting the eruptive fissure to the deep dike. A similarly narrow path is also postulated for other volcanoes such as Piton de la Fournaise [Fukushima et al., 2010; Traversa et al., 2010].

6.5. Eruption Model

[60] In order to distinguish whether the deep or the shallow dike has the most influence on the opening of the other one, we use the MBEM to compute normal stress changes induced by each inverted dike on the other. For these computations, we use the uncorrected overpressures P₀ = 0.9 MPa and P₁ = 0.4 MPa determined for the shallow and deep dike (Tables 3a and 3b), respectively, considering a homogeneous medium with a Young's modulus E = 5 GPa, as in the inversions. Induced normal stresses are computed at the centroid of each triangular boundary element. Figure 11 shows that each dike decreases the normal stresses of the other. The shallow dike mainly decreases the normal stresses on the upper, northern part of the deep dike (Figure 11a), while the deep dike mainly decreases normal stresses on the southern part of the shallow dike in an area extending from the bottom of the dike to the surface (Figure 11b).

[61] Geophysical observations [Battaglia et al., 2005], numerical [Dahm, 2000] and analog models [Rivalta and Dahm, 2005] indicate that dikes grow upward or laterally from the source region. Thus it is unlikely that the stress changes induced by the shallow dike would promote the injection of magma downward to the deeper dike. Conversely, the stress decrease induced by the deep dike could favor the injection of magma from the deep dike up into the shallow dike. The area of maximum stress decrease is coincident with the area where radioactive disequilibria indicate a different magmatic source to that of the northern part of the eruptive dike [Tedesco et al., 2007]. Moreover, the deep dike induces a 0.07 MPa pressure decrease in the host rocks at the level of the shallow reservoir, corresponding to an overpressure increase in the reservoir. If the reservoir is at a critical equilibrium and connected to the summit lava lake, as indicated by the petrographical data [Demant et al., 1994], this stress change might be sufficient to initiate the reservoir and magma column failure, triggering the intrusion of a shallow dike.

[62] A likely scenario (Figure 12) is that the eruption was initiated by the supply of magma to a deep reservoir, ten months before the eruption, as suggested by the increased number of long period earthquakes and the onset of continuous harmonic tremor before the eruption [Kavotha et al., 2002/2003]. Such an early supply is also consistent with the model of lava flow rheology which indicates that lavas emitted on the southernmost fissure had probably been cooled in the weeks or months before the eruption [Giordano et al., 2007]. This deep reservoir might have subsequently fed the deep dike.

[63] As indicated by our stress analysis, the intrusion of the deep dike could have promoted tensional failure of the shallow reservoir. As this reservoir is permanently connected to the lava lake [Demant et al., 1994], failure of the reservoir could have induced the opening of the shallow dike and the subsequent draining of the lava lake.

[64] Lavas from the lava lake were emitted first as they were closer to the surface and their magmastatic pressures promoted dike opening, followed by lavas from the shallow reservoir. Draining of the shallow reservoir could have induced failure of the top of this reservoir along ring dikes, leading to the unusual earthquakes recorded during the eruption [Shuler and Ekström, 2009]. Even if the deep intrusion was initiated first, it might only have breached the surface after the shallower intrusion since it was further from the surface. Since fracture toughness is depth-dependent [Rivalta and Dahm, 2006] the deeper intrusion probably propagated to the surface more slowly than the shallow intrusion. The presence of the deep dike could have facilitated the connection between the deep dike and the southern part of the shallow dike, leading to the eruption of magma from the deep source south of the eruptive fissure. This scenario involving the lava lake, a shallow reservoir and a deep reservoir for the southernmost fissure is consistent with the geochemical analysis and model of Tedesco et al. [2007].

[65] Injection of magma into a deep dike could have triggered the precursory and post-eruptive short- period earthquakes [Kavotha et al., 2002/2003; Tedesco et al., 2007] by inducing slip on faults close to the dike tip which were on the verge of failure [Rubin, 1992]. This has often been observed during dike intrusions, for instance in Afar in 2006 [Grandin et al., 2009; Keir et al., 2009], at the volcanic island of Izu in 2000 [Nishimura et al., 2001] or during the 2007 Natron faulting and diking event [Baer et al., 2008; Calais et al., 2008; Biggs et al., 2009]. Large earthquakes, such as the 3 February 2002 5.9 M_w Bukavu event [d'Oreye et al., 2011], or the October 24, 2002 6.2 M_w south Lake Kivu event (Figure S8 in the auxiliary material), might remain undetected by interferometry if they were to occur beneath Lake Kivu or in areas of low coherence.

[66] On the 9 descending and 8 ascending available RADARSAT-1 interferograms that cover the year before the eruption, we do not see any evidence of pre-eruptive inflation of a deep reservoir, but it might go undetected if it was too deep or associated with only a small volume increase, or because the area affected by the magma injection was incoherent, as would be the case if it was located beneath Lake Kivu.

[67] The lack of seismic and continuous geodetic measurements prevents us from tracking the progression of magma at depth. Therefore, compared to better-studied rifting events, such as the Manda-Hararo 2005–2010 diking events [Keir et al., 2009; Grandin et al., 2010] or the Puu' Oo 1983 eruption [Rubin et al., 1998], it is not possible to detect if the deep dike was fed by a mid-segment reservoir.

6.6. Comparisons With the 1977 Eruption

[68] The 2002 eruption has many common characteristics with the 1977 event. Both eruptions are associated with rift-parallel fissure networks, indicating that the intrusion directions are related to the direction of rift extension.

[69] The two estimated erupted lava volumes are similar (Table 4, column 4), and, in both eruptions, the summit lava lake drained. In 1977 the lava lake level was about 155 m below the crater rim, whereas it was about 250 m below the rim before the 2002 eruption (Table 4, column 1). Not only was the lava lake level higher in 1977 than in 2002, but the diameter of the drained crater was wider in the former [Durieux, 2002/2003], leading to a much larger estimate for the drained volume in 1977 (234 × 10⁶ m³) [Tazieff, 1977] compared to 2002 (60 × 10⁶ m³) [Komorowski et al., 2002/2003] (Table 4, column 2).

[70] For the 1977 eruption, an intruded volume of 212 × 10⁶ m³ was estimated from the difference between the drained lake and the lava flow volumes [Komorowski et al., 2002/2003; Tazieff, 1977]. Following the same approach, an intruded volume of 60 × 10⁶ m³ was estimated from the drained volume in 2002 (Table 4, column 3) [Komorowski et al., 2002/2003]. However, inverting for InSAR data, we found that both a shallow and a deep dike were involved in this eruption, leading to an intruded volume of 210 × 10⁶ m³, more than three times larger than this estimate. Our modeling is consistent with the lava lake draining into the shallow dike. This shows that, when estimating magma volumes, volumes drained from the lava lake should be considered as minimum volume estimates as they do not account for magma supplied by reservoirs.

[71] The 1977 eruption appears to have induced a 2.4-fold increase in eruption rates at neighboring Nyamulagira volcano [Burt et al., 1994; Wadge and Burt, 2011], while the 2002 eruption may have decreased the lava output rate at Nyamulagira [Wadge and Burt, 2011]. Large dikes, such as the ones in this study, imply stress perturbations at a rift scale that could also affect the 15 km-distant Nyamulagira volcano and change its eruption rate or its intruded versus extruded magma volume. Deformation associated with the 1977 eruption was not measured, so the exact perturbation related to this event cannot be investigated and compared with that of the 2002 eruption.

6.7. Isotropic Lithostatic Crustal Stresses

[72] The deep dike overpressure is smaller than values estimated in areas where extension is accommodated by normal faulting. Horizontal stresses needed for slip on faults correspond to a critically stressed crust with hydrostatic pore pressure [Rubin, 1990; Buck, 2006; Townend and Zoback, 2000]. In this case, Rubin [1990] estimated that an average overpressure of 54 MPa could be reached at 5 km depth. However, the overpressures we determined are close to values estimated when the stress state in the crust is lithostatic and isotropic, such that the horizontal compressive stress equals the rock pressure. If the rock stress close to the dike was lithostatic and isotropic, the overpressure at the center of the dike z_d (corresponding to the mean overpressure), would be the difference between the magma and rock pressure:

where z_max and z_min are the maximum and minimum depths of the dike, ρ_m is the magma density and ρ_r(z′) is the rock density, linearly varying with depth, as indicated in Figure 6b. g is the acceleration of gravity.

[73] For the estimation of the magma density, we consider that the magma is degassed. The InSAR displacements cover a six-month period extending from the eruption until mid-July 2002, during which time the SO₂ emission rates were low. This led Carn [2002/2003] to assume that the emitted magma had been degassed during the long residence time in the crater and upper conduit of the volcano. For the degassed magma density in the shallow dike, we take a value of ρ_m = 2650 kg/m³, which is the nephelinite density at Nyiragongo [Carn, 2002/2003], and for the magma density of the deep dike we take ρ_m = 2800 kg/m³.

[74] Equation (8) allows us to determine that the overpressure corresponding to lithostatic isotropic stresses are 1.1 MPa, and 4.5 MPa for the shallow and deep dikes, respectively. These values are close to the overpressures estimated from the inversion of ground displacements, which are 0.9 MPa for the shallow dike (Table 3a) and 7.4 MPa for the deep dike, after correcting for the low Young's modulus used for the inversions (15). Our analysis indicates that, close to the dikes, stresses are incompatible with the tensional tectonic stresses needed to produce stretching via normal faulting [Buck, 2006; Rubin, 1990; Townend and Zoback, 2000].

[75] If magma pressure was higher than the hydrostatic pressure assumed in equation (8), then the overpressure would be increased. The low overpressure we determined thus indicates that the magma pressure is low. This result is consistent with the observation that, in some places, magma was visible in the eruptive fissure without being actively erupted.

[76] Few studies of diking events document the overpressure, as most of them rely on kinematic inversions of geodetic data. Stress changes were determined for the different episodes of the 2005–2009 Afar rifting cycle by Grandin et al. [2010]. These authors did not invert surface displacements for stress changes but computed stress changes a-posteriori from the opening distribution, which is such that it represented the best compromise between the data fit and the smoothness of the dike opening. No assumption was made about the uniformity or depth dependence of overpressures, as would be expected for magma-filled fractures. Obtained stresses were very large at the center of the open patches, reaching 30 MPa for a 5 km-deep dike. We attribute the large overpressures obtained by Grandin et al. [2010] to the method they used for stress determination.

6.8. Influence of the Rift Extension on the Dike Directions

[77] Both dikes are subvertical and more or less perpendicular to the rift extension direction (Figure 1) [Stamps et al., 2008], the shallowest dike deviating more from the extension direction than the deeper dike. Northwest of the Nyiragongo edifice, the shallower dike follows the N20°W weakness direction [Meyer, 1953] linking Nyamulagira and Nyiragongo, while further south the dike follows the north-south direction taken during the course of the 1977 eruption. This direction is intermediate between the N20°W weakness direction and the N20°E rift direction, suggesting the combined effect of the two on the eruptive fissure emplacement. The fact that the direction of the deeper dike is perpendicular to the extension direction indicates that this dike direction is guided by the rift extension. As crustal stresses are found to be isotropic, it is likely that this direction is favored because of a reduced tensile strength inherited from previous dike intrusions rather than because it is a direction of minimum principal stress.

[78] The deeper dike is aligned to several other rift-related structures. For instance, it is aligned to the Idjiwi horst (Figure 2a), located south of Lake Kivu. Wong and Von Herzen [1974] showed that the horst structure extends up to the northern shore of the lake. The bounding fault of the horst could have guided the magma path. The deep dike could be connected to an undated Ngangi eruptive fissure (Figure 2b), 10 km west of the 2002 fissure, which strikes N8°E. This fissure is associated with the alignment of several eruptive cones and hornitos, and has localized uplift of about 10 m, as confirmed by a field survey in January 2010. At the time of the January 2002 eruption, the local population reported smoke emissions along part of the Ngangi fissure. Eruptions beneath Lake Kivu are probably common since bathymetric data [Capart, 1955; Schmid et al., 2010] indicate the presence of numerous old volcanic cones on the lake floor (Figure 2a). The southernmost cones are aligned with the deep dike, as determined from our inversions, indicating that this structure might have been active in the past and breached the lake floor.