Magma Pathways and Their Interactions Inferred from InSAR and Stress Modeling at Nyamulagira Volcano , D . R . Congo

A summit and upper flank eruption occurred at Nyamulagira volcano, Democratic Republic of Congo, from 2–27 January 2010. Eruptions at Nyamulagira during 1996–2010 occurred from eruptive fissures on the upper flanks or within the summit caldera and were distributed along the ~N155E rift zone, whereas the 2011–2012 eruption occurred ~12 km ENE of the summit. 3D numerical modeling of Interferometric Synthetic Aperture Radar (InSAR) geodetic measurements of the co-eruptive deformation in 2010 reveals that magma stored in a shallow (~3.5 km below the summit) reservoir intruded as two subvertical dikes beneath the summit and southeastern flank of the volcano. The northern dike is connected to an ~N45E-trending intra-caldera eruptive fissure, extending OPEN ACCESS Remote Sens. 2015, 7 15180 to an ~2.5 km maximum depth. The southern dike is connected to an ~N175E-trending flank fissure extending to the depth of the inferred reservoir at ~3.5 km. The inferred reservoir location is coincident with the reservoir that was active during previous eruptions in 1938–1940 and 2006. The volumetric ratio of total emitted magma (intruded in dikes + erupted) to the contraction of the reservoir (rv) is 9.3, consistent with pressure recovery by gas exsolution in the small, shallow modeled magma reservoir. We derive a modified analytical expression for rv, accounting for changes in reservoir volume induced by gas exsolution, as well as eruptive volume. By using the precise magma composition, we estimate a magma compressibility of 1.9–3.2 × 10 Pa and rv of 6.5–10.1. From a normal-stress change analysis, we infer that intrusions in 2010 could have encouraged the ascent of magma from a deeper reservoir along an ~N45E orientation, corresponding to the strike of the rift transfer zone structures and possibly resulting in the 2011–2012 intrusion. The intrusion of magma to greater distances from the summit may be enhanced along the N45E orientation, as it is more favorable to the regional rift extension (compared to the local volcanic rift zone, trending N155E). Repeated dike intrusions beneath Nyamulagira’s SSE flank may encourage intrusions beneath the nearby Nyiragongo volcano.


Introduction
It is a challenge to quantitatively predict the position and direction of dike intrusions and resulting eruptive fissures at volcanoes, because they are governed by the interplay between several factors, such as a heterogeneous regional stress field, preexisting discontinuities and heterogeneous and anisotropic properties of rocks.Nevertheless, dikes at basaltic volcanoes usually show some regularity in their emplacement pattern, as eruptive fissures and vents tend to align with the direction of the maximum regional principal stress [1].For instance, radial and circumferential diking alternate at the Galapagos volcanoes [2][3][4]; persistent diking occurs along well-developed volcanic rift zones at Kilauea Volcano (e.g., [5]); and summit eruptions alternate with proximal and distal eruptions at Piton de la Fournaise volcano [6].
The 2010 Nyamulagira eruption was observed by Interferometric Synthetic Aperture Radar (InSAR) with different acquisition geometries, allowing detailed numerical modeling of the involved sources of deformation.Although the event was a rather typical upper flank eruption, it may have perturbed Nyamulagira's plumbing system significantly, as the following eruption in 2011-2012 occurred ~12 km ENE of the summit (Figure 1) [10].The mechanism responsible for such a change in eruption location and behavior is unclear.
Figure 1.Geological setting of the Virunga Volcanic Province (VVP), western branch of the East African Rift System (EARS), north of Lake Kivu.The rift extension direction in this area is ~N110E, with an estimated extension rate of 2-2.8 mm/year [11,12].The Goma-Gisenyi urban area is outlined in white on the northern shores of Lake Kivu.Eruptive fissures and lava flows (see the legend) were mapped from radar and optical images [7].Nyamulagira's caldera rim is outlined with a white dashed line.Faults are from [13], and rose diagrams of VVP's faults (yellow) and eruptive fissures (red) are from [9].The white rectangle gives the extent of Figure 2.
In this study, we present co-eruptive InSAR observations of the 2010 eruption at Nyamulagira volcano.We invert the InSAR deformation field to infer volume changes, geometries and locations of sources of deformation involved in the eruption.Finally, we examine the normal stress changes induced by the modeled 2010 sources of deformation on Nyamulagira's rift zones and provide insights about the relationships between recent eruptions and changes in eruption behavior.
Figure 2. Data, 1st column: most coherent wrapped interferograms covering the entire 2010 Nyamulagira eruption (heights of ambiguity, i.e., the topographic altitude difference, which generates one residual topographic fringe in a differential interferogram [14]; for interferograms, (a)-(d) are 45, 430, 704 and 3000 m, respectively).Pixels with coherence lower than 0.2 have been masked out (the unmasked wrapped interferograms are shown in Figure S1, and the four selected unwrapped interferograms are shown in Figure S2).Nyamulagira's summit caldera is outlined with a dashed white line.Note that the ALOS interferogram has been rewrapped to the ENVISAT's wavelength of 5.6 cm.Model, 2nd column, and residuals, 3rd column, for the preferred best-fit model (Figure 3), including two dike intrusions associated with the 2010 eruptive fissures (bold green lines) and a deflating reservoir, for the four interferograms (a)-(d).

Geologic Setting
Nyamulagira is a basaltic shield volcano located in the Virunga Volcanic Province (VVP), Democratic Republic of Congo.The VVP volcanism is probably controlled by the transfer zone that links the Lake Edward and Lake Kivu segmented extensional basins located along the western branch axis of the EARS [15][16][17][18].Nyamulagira is one of Africa's most active volcanoes, with eruptive episodes occurring every 1-4 years since the 1980s [10].Eruptive centers and fissure distributions mostly follow two orientations: ~N155E and ~N45E [9].The ~N45E orientation corresponds to the strike of rift-related structures in the VVP (see the fault and eruptive fissure rose diagrams in Figure 1).The ~N155E orientation crosses the caldera, extends towards Nyiragongo's edifice and can be considered as a volcanic rift zone.

Characteristics of Recent Eruptions
Eruptions at Nyamulagira during 1996-2010 occurred from eruptive fissures on the upper flanks or within the summit caldera and were distributed along the ~N155E rift zone (Figure 4).The 1996The , 1998The , 2000The , 2004The , 2006 and 2010 eruptions correspond to typical upper flank (<3 km from the caldera) short-lived eruptions, which last a few days to weeks and generate between ~45 and 70 × 10 6 m 3 of lava [9].The magma that feeds such upper flank eruptions likely intrudes as subvertical dikes from a storage area located at a depth of up to 5 km below the ground surface [19,20].In contrast, the 1989, 1991-1993, 2001 and 2011-2012 eruptions correspond to less frequent (roughly decennial) long-lived eruptions, which last from several months to years, produce larger volumes of lava (i.e., >80 × 10 6 m 3 ) and occur in the summit area (as in 2001) or >9 km from the summit [9].The location of the ~N70E 2011-2012 eruptive fissures (Figure 4) corresponds to the eastern extension of an unusual ~4 km-deep long-period (LP) earthquake swarm oriented ~N45E, which was recorded the month preceding the 2010 eruption (Figure 5).Dikes feeding these larger volume eruptions might initiate from a deeper (~20-km depth) magma reservoir located between Nyamulagira and Nyiragongo [9].An alternative hypothesis is that eruptions fed by dikes parallel to local EARS structures (i.e., oriented ~N45E, Figure 1) have a longer duration and larger volumes because they are more favorably oriented with respect to the remote stress associated with the rift extension [8].   2. See [10] for a detailed review of the four eruptive stages.SP and LP stand for short and long period, respectively.

The January 2010 Eruption
On the morning of 2 January 2010 at 2:15 a.m.(local time, UTC/GMT+2), a swarm of short-period volcano-tectonic earthquakes and associated volcanic tremors began at Nyamulagira and is interpreted to indicate the eruption onset [10].Clear seismic precursors (e.g., increased number of long-and short-period events, along with increased sequences of volcanic tremors) started only 2 h before the onset of the eruption [10].Time series of deformation computed from InSAR data [21] indicate pre-eruptive vertical (uplift) displacements at the summit (5 cm) and close to the 2010 eruptive fissures (3 cm).After the start of the eruption, seismicity was mainly characterized by continuous volcanic tremors (Figure 5).Intense lava fountains were observed within Nyamulagira caldera, in a pit crater and along ~N45E-trending fissures (Figure 4).Furthermore, lava fountains were also observed along a new ~600 m-long ~N175E fissure on the southern flank of the volcano, ~1.5 km from the caldera rim, building a new volcanic cone.After ~8 days, the eruptive activity was localized to this circular vent on the southern flank (Figure 5; see [10] for a detailed review of the four eruptive stages).The main lava flow reached up to ~10.6 km from the eruptive vent, covering portions of the 2006 lava flows (Figure 1) before the eruption ceased on 27 January 2010.Due to the poor configuration of the seismic network around the volcano and the use of an over-simplistic velocity model [10,22], the location of earthquakes is poorly constrained (horizontal location uncertainties are in the order of ~10-15 km, and the vertical uncertainty is likely even larger).Therefore, we only consider swarms and general trends to be significant [10].

InSAR Processing
Data from the European Space Agency (ESA)'s ENVISAT satellite and from the Japanese Exploration Agency's ALOS-1 satellite imaged the deformation of Nyamulagira during the January 2010 eruption.The most coherent interferograms covering the entire eruption with the shortest time span possible are shown in Figure 2 (Figure S1 in the Supplementary Material shows all of the unmasked interferograms available).The interferograms were processed using the JPL/Caltech ROI_PAC SAR software, multilooked with a factor of 2 in range and 10 in azimuth and filtered with an adaptive filter to optimize the signal-to-noise ratio [23].We removed the orbital and topographic interferometric phase contribution using precise orbits and a 30-m-resolution digital elevation model generated by the NASA Shuttle Radar Topography Mission [24].The interferometric phase was unwrapped (Figure S2) with the SNAPHU algorithm [25].Line of sight (LOS) ground displacements of up to ~15 cm are visible northwest of the caldera eruptive fissure and correspond to a range decrease (ground displacement toward the satellite) in the ascending interferograms (Figure 2a,b and  d) and a range increase (ground displacement away from satellite) in the descending interferogram (Figure 2c).A range decrease of ~6 cm is visible north of the caldera fissure in the descending interferogram only.A range decrease of ~25 cm is visible west of the southern flank fissure in the ALOS ascending interferogram (note that this area is incoherent in the other interferograms).Deformation up to ~12 cm is visible east of the southern flank fissure and corresponds to a range decrease in the ascending interferograms and a range increase in the descending interferogram.Finally, a few centimeters of range increase is also visible in the descending interferogram north of the southern flank fissure.Note that the maximum amplitude of horizontal deformation (<25 cm, [21]) is not large enough to be resolved by offset-tracking methods, for which the errors on the offset estimations will be of the same order as the deformation (e.g., [26,27]).Moreover, the deformation field is dominated by the ~EW displacements [28] induced by the dike intrusions (Figure S2), making an azimuthal offset calculation unlikely to constrain our inversions and models further.

Boundary Elements Modeling and Non-Linear Inversions
We modeled the deformation for the 2010 Nyamulagira eruption using the 3D-mixed boundary element method [29].This numerical modeling approach explicitly takes into account the topography, mechanical interactions between sources and allows realistic dike models connected to the eruptive fissures [30].The convention for stresses follows rock mechanics' convention where compression is positive.We assume a Poisson's ratio of 0.25.The Young's modulus is assumed to be 5 GPa, corresponding to the value inferred from seismic velocities for the upper few kilometers of the crust in the VVP [31].
Prescribed boundary conditions are stress vectors.They are null at the ground surface.Sources are assumed to undergo constant pressure changes, corresponding to the difference between the magma pressure and the normal stress exerted by the host rock.Although the way to include gravitational stresses and other remote stresses is a matter of debate [32,33], we chose to neglect vertical gradients corresponding to gravitational forces, as a previous study conducted for dikes at Piton de la Fournaise volcano [30] showed that vertical gradients could not be resolved.For dikes, the constant overpressure boundary condition is consistent with hydraulic connection within all parts of the dikes.Dikes are parameterized with 6 geometrical parameters allowing for physically-meaningful parameter ranges (see Supplementary Methods and Figure S4 in the Supplementary Material), and reservoirs are parameterized with 4-8 parameters, depending on the tested geometry.These assumptions are consistent with [10], who inferred that the two 2010 eruptive fissures opened at the same time and were fed from the same magma ascent.
Ground surface element size affects ground surface displacement patterns, while source element size affects ground surface displacement amplitudes, and thus, the stress changes estimation [30].Consequently, two meshing strategies are followed: one for the ground surface and one for the source meshes.For the ground surface, mesh elements have sizes that increase away from the maximum of displacement (Figure S5 [34]).The element size, which provides the best compromise between accuracy and computation time, is determined by successively decreasing the elements sizes until surface displacement is unchanged.To limit computation time, for deformation sources, coarse elements are used in inversions, and after inversions have converged, final computations with 4-times more elements are performed to adjust the source pressure changes.
As surface displacements are a non-linear function of the sources' geometry and location, nonlinear inversions of the deformation sources are required.To invert the InSAR data and solve for the most likely models, we used a nearest-neighbor algorithm [35].To make the inversions computationally manageable, we subsample the unwrapped interferograms following a circular subsampling method [30,36], leading to a total of 2118 subsampled points.
Generally, when inverting for model parameters, increasing the number of model parameters decreases the misfit; however, the lowest misfit model determined is not necessarily the most likely.In order to select the most likely model, we select the model that minimizes the Akaike information criterion (AIC) [37,38], which is defined as: where  is the number of inverted parameters, N is the number of subsampled data points and  2 is the misfit.

Numerical 3D Modeling Results
The January 2010 co-eruptive InSAR deformation is the result of successive eruptive processes that occurred between mid-December 2009 and early February 2010.The InSAR data provide a "snapshot", which includes several potential volcanic processes: a pre-eruptive reservoir inflation, the intrusion of magma in dikes and sills, a co-eruptive reservoir deflation and a post-eruptive pressure restoration from gas exsolution (Figure 5).Unfortunately, the lack of continuous ground-based measurements means that these potential events cannot be observed individually, so their occurrence is a matter of debate and will be discussed in Section 5.
To model the January 2010 co-eruptive deformation, we first inverted the InSAR data for a single dike connecting the two eruptive fissures and found that it could explain only 68% of the observed deformation (Table 1).Given the localized patterns of deformation observed (Figure 2), the January 2010 eruption must have been fed by two dikes, spatially close to each other.Following [39], we assume that they interact, and we thus invert them simultaneously.Considering two distinct dikes, associated respectively with the summit and the southern flank eruptive fissures, only slightly improves the fit with 71% of the observed deformation explained at best (Table 1).Moreover, adding a simple quadrangular connection between both dikes beneath the southern flank provides the worst fit, with only 60% of the observed deformation explained at best (Table 1).Even with two dikes, the circular deformation signal north of the caldera, corresponding to a range increase visible in the ENVISAT descending interferogram (Figure 2c), remains unexplained.Such a signal may however be caused by the pressure change in a reservoir located beneath the caldera.Hence, we then invert for one or two dike(s) plus one of the following sources: a circular sill, a rectangular sill, a spherical reservoir or a laccolith dome-shaped intrusion (Table 1), assumed to experience a uniform static pressure change.When modeling magma reservoirs, there are trade-offs between the reservoir dimension and the pressure change, and only the latter can be determined with confidence.Additional data, such as magma extrusion rates, are required to resolve this trade-off [40].Here, we arbitrarily fixed the reservoir radius to 500 m.In Section 5.2, we will discuss the consistency of this reservoir dimension with the erupted volume and host rock failure ability.
The best-fit model, giving the lowest RMS error and AIC value and explaining 90% of the observed deformation (Table 1, Figure 2), includes two dikes and a spherical deflating magma reservoir located beneath the caldera.With this model, the dike connected to the northern eruptive fissure (the "northern dike") in the caldera is subvertical, extends to ~1.5 km beneath the summit and has a horizontal bottom side (Figure 3 and Table 2).It experiences a volume increase of ~0.74 × 10 6 m 3 .The second dike, connected to the southern eruptive fissure on the SE flank (the "southern dike"), extends to a greater depth, reaching ~3.5 km beneath the surface.It has an ~6 km-long bottom side dipping at ~40 degrees northward (Figure 3 and Table 2).It experiences a larger volume increase of ~4.2 × 10 6 m 3 .The spherical magma reservoir, located beneath the SW depression in Nyamulagira caldera (Figure 3) and centered at a depth of ~4 km beneath the summit, experiences a volume decrease of ~5.1 × 10 6 m 3 .

Table 2.
Best-fit parameters for the most likely model using two dikes and a spherical reservoir.The 95% confidence intervals on the 12 inverted parameters are shown (Figure S6, [41]).

Parameter
Dike Connected to the Caldera Fissure Dike Connected to Southern Flank Fissure

Sources of Deformation
The preferred best-fit model, including two dikes and a spherical reservoir, leads to reasonable residuals (Figure 2).There are some residuals in the summit area, probably advocating for a more complex magma reservoir/northern dike geometry or a more complex rheology.However, it is difficult to further constrain the geometry of the summit deformation source due to the lack of coherence on the vegetated flanks.Similarly, some residuals are visible east of the southern flank eruptive fissure, possibly because of new ash deposits.The bottom of the southern dike extends beneath the southern flank of Nyamulagira where the opening and subsidence of a graben structure is visible [20], as commonly observed at the upper tip of blade-shaped dikes [42][43][44][45].Therefore, there is probably a pre-existing area of weakness close to the surface in which magma intrudes laterally along dikes to feed the SE flank eruptions from the shallow summit reservoir beneath the caldera.The upper southeastern flank of Nyamulagira may be developing a volcanic rift zone, similar to the well-developed volcanic rift zones observed on the volcanoes of Hawaii, La Ré union and the Canary Islands (e.g., [46]).
As both fissures were fed by the same batch of magma [10], the dikes and reservoir are probably connected, but our modeling method cannot resolve this connection.The connection between both dikes is probably via another dike, as this is the most efficient means of transporting magma through the lithosphere [47,48].Such a dike may be too narrow to cause detectable InSAR displacements.Indeed, the data fit breaks down when we try to add a sizeable connecting dike (Table 1), but if the connecting dike has a width smaller than 400 m at a 1-km depth, forward 3D-MBEM models show that it is undetectable.This dike width is similar to the eruptive fissure width, indicating that flow in the connecting dike could be the same as the flow at the surface.A similarly narrow path is also postulated for other volcanoes, such as Nyiragongo [31] and Piton de la Fournaise [39,49].
The fissure that opened on the caldera floor follows an ~N45E orientation, which contrasts notably with the N155E orientation of the main rift zone and the southern flank eruptive fissure.However, this N45E orientation corresponds to a fracture system that opened during the 1938-1940 major eruption [10].Furthermore, during that same eruption, the SW caldera depression formed as the result of a collapse of the caldera.The magma reservoir modeled in this study is located in the same zone as the reservoir implied in the 1938-1940 eruption [10].In 2010, this reservoir could have induced the same stress field as in 1938-1940, reactivating N45E-trending pre-existing structures.The eruptive fissures are not located at the crest of the reservoir, as expected when a spherical reservoir fails under tension [32], but away from the summit.This might be explained by host rock stresses inherited from previous eruptions.

Magma Reservoir Size
Assuming the reservoir is spherical, the maximum overpressure sustainable before tensile failure of the host rock is given as a function of the lithostatic stress, Pl, and the tensile strength, T0 [39,56], as: where   =    , with   indicating host rock density,  indicating the acceleration of gravity and  indicating the depth to the top of the reservoir.If we assume   = 2600 kg/m 3 [57],  0 = 0.5-5 MPa [58] and use a value of  = 3.5 km from our inversions, the overpressure ∆  = 179-183 MPa.If we consider that the pre-eruption excess volume   , which induced the failure, corresponds to the pre-eruptive inflation detected by InSAR,   =5 cm [21], we can estimate the radius of this reservoir.
Assuming the reservoir volume change corresponds to that of a spherical reservoir in an infinite medium (the reservoir depth/radius ratio is >5), we can compute this radius from [59]: where  and  are Young's modulus and Poisson's ratio, respectively, and   =    2 (1 − ) ⁄ .Because of the 1/3 power, the estimated volume has little sensitivity to the overpressure or the volume change.We find that the reservoir radius should be smaller than R = 250 m in order for failure to occur.Since InSAR measurements are not continuous, this reservoir radius is likely larger.For instance, if pre-eruptive inflation corresponds to the co-eruptive deflation, then the reservoir should be smaller than R = 315 m.These values are of the same order as the assumed reservoir radius, suggesting that our model is plausible.As surface deformation measurements are not sensitive to the reservoir radius for depth/radius ratios > 5, a change in the reservoir radius will not affect the depth and the reservoir volume changes determined from the inversions.
No reactivation of the ring fault around the summit caldera was observed.This occurs when magma reservoirs undergo under-pressure.Here, the depth of the eruptive fissure is at a higher elevation than the reservoir, and at the end of the eruption, degassed lava was emitted, implying that the reservoir is probably at lithostatic equilibrium at the end of the eruption.

Magma Budget
The InSAR data covering the Nyamulagira 2010 eruption require a deflating source beneath the caldera, which we modeled as a spherical reservoir at ~3.5 km beneath the summit (Figure 3) experiencing a volume decrease of ~5.1 × 10 6 m 3 .The intruded magma volume in the two dikes is ~4.9 × 10 6 m 3 , similar to the reservoir volume decrease.The estimated extruded lava flow and tephra volume is 45.5 ± 7.6 × 10 6 m 3 [10] and 10 × 10 6 m 3 [8], respectively, leading to an average dense rock equivalent volume of ~42.3 × 10 6 m 3 (assuming a lava, tephra and dense rock density of 2200, 1000 and 2600 kg/m 3 [8], respectively).The total emitted volume is thus ~47.2 × 10 6 m 3 (= the sum of dike volumes + dense rock equivalent erupted volume), leading to a ratio (rv) between the total emitted volume and the reservoir volume decrease of ~9.3.This ratio is an intermediate value between commonly-reported values in the EARS, for instance of ~20-40 [45] for the Dabbahu intrusion and of ~3 for the Dallol intrusion [60].Refill of the shallow reservoir from a deep reservoir is not necessarily required to explain this discrepancy, as the pressure in the reservoir could have been restored by gas exsolution in the magma chamber [50].
We follow an approach similar to the one of Rivalta and Segall [50] to express rv.We assume that mass is conserved between magma emitted from the reservoir and magma transmitted to the edifice and that magma is compressible.However, we additionally account for the lava erupted at the surface (see Appendix A) and assume that pressure in the magma chamber before the eruption,   , corresponds to the maximum sustainable overpressure defined above,   =   +   = 3  +   , and that, after the eruption, pressure in the magma chamber returns to a lithostatic state.Thus, we obtain the following: where   and   are the magma and reservoir compressibility,   is the dense rock equivalent volume of emitted products,   is the reservoir volume change,   is the dense rock equivalent density of emitted products (lava flows + tephra) and   (  ) and   (  ) are the densities of magma in the reservoir after and before the eruption.Reservoir compressibility can be computed by assuming a spherical reservoir and neglecting the influence of the ground surface [50], leading to   = 3 (1 + ) (2) ⁄ .Compressibility of magma,   , is derived from   = 1   ⁄     ⁄ , assuming   remains constant over the pressure interval, so that: The variation of magma density with pressure depends on the amount of exsolved gas, given by: where   () is the mass fraction of gas in the magma and   ()and   () are the densities of gas and liquid, respectively, all of which vary with pressure.  () can be determined from the analytic expression derived from Henry's law [61] or derived from more complex models [62][63][64].The gas density is given by the ideal gas law   () =    ⁄ , where M is the gas molar mass, R is the gas constant and T is the absolute temperature, while the liquid density takes the liquid compressibility,   = 1 − 2 × 10 −10  −1 , into account:   =   (1 +   ), where   = 2500 kg m 3 ⁄ is the density of the liquid phase of magma at atmospheric pressure.
Here, the volatiles considered are CO2 and H2O, and the magma compositions are estimated from recent studies of the volcano [10,65] (Table 3).Using the online code for thermodynamics computation of Papale et al. [64], we find that for a 1200 °C magma at a 3.5-km depth, almost all CO2 is gaseous and all H2O is liquid, so that   2 () = 0.1 wt% and   2  () = 0 wt%.At the pressure corresponding to rupture of the reservoir,   , we find   2 (  ) = 0.05 wt%.Successively using Equations ( 5) and ( 6), we get a magma compressibility of   = 1.9-3.2× 10 9 Pa −1 .If we consider a reservoir compressibility for a source in an infinite medium [50] and take   = 2600 kg/m 3 for the dense rock equivalent density and   = 47.2 × 10 6 m 3 , using Equation ( 4), we obtain rv = 6.5-10.1, a value consistent with that determined from modeled volume changes.The contribution of the emitted lava flow to rv is 0.5, thus it accounts for approximately only 5% of rv.This value of   indicates that the magma emitted can be entirely issued from the 3.5-km depth shallow magma chamber and that discrepancy between the emitted and the reservoir volume change are compensated by CO2 exsolution.

Stress Analysis
To understand the possible origin of the 2011-2012 eruption, we compute static stress changes induced by the modeled 2010 eruption deformation sources.Static stress changes can be used as an indicator for dike unclamping (e.g., [31,[66][67][68][69][70]). Indeed, each dike intrusion changes the normal stress on pathways for potential diking.The sign of the induced normal stress change indicates whether subsequent intrusions are favored or prevented.If the normal stress changes are negative, subsequent dike emplacement is favored (unclamping); otherwise, subsequent dike emplacement is prevented (clamping).Note that we do not intend to forecast the shallow feeder dike orientation for the 2011-2012 eruption, but are interested in understanding the initiation of dikes.To quantitatively estimate the stress changes induced by the 2010 Nyamulagira eruption, we compute normal stress changes on five potential planes induced by the best-fitting deformation sources.From the distribution of eruptive vents, fissures and pre-existing rift-related structures (Figure 1), five potential vertical dike surfaces were considered (Figure 4): an ~N155E direction corresponding to the more active rift zone, an ~N45E direction crossing the summit caldera and an ~N45E south caldera (SC) direction, as well as an ~EW direction and an ~N65E direction connecting to the 2011-2012 eruptive fissure corresponding to a radial dike emplacement.Potential dike surfaces investigated range from the ground surface down to three kilometers below sea level (b.s.l.).
The 2010 eruption deformation sources induce unclamping in all potential dike surfaces (Figure 6 and Figure S3).Along the ~N45E surface crossing the caldera, the ~EW and ~N65E radial investigated surfaces, the unclamping is only significant in the shallowest part of the surfaces (<~2 km beneath the ground surface) and in the immediate proximity of the modeled reservoir.The ~N155E and the ~N45E SC surfaces are unclamped over large areas (66 and 111 km 2 , respectively) to greater depths.Stress changes from the 2010 intrusions and reservoir beneath the SSE flank could have thus encouraged the vertical ascent of magma from a deeper crustal reservoir along the ~N45E SC surface below the volcanic edifice.This interpretation is consistent with the hypothesis of a deep (~20-30-km depth) magma reservoir feeding the less common distal lower flanks' eruptions [19].We also calculated the normal stress change induced by the 2010 Nyamulagira eruption on the ~NS potential dike surface at Nyiragongo, determined from the modeling of the 2002 Nyiragongo eruption [31] (Figure S7).The normal stress change is slightly negative close to Nyiragongo summit (the minimum value is −0.0015MPa).Thus, dike intrusions in the SSE flank of Nyamulagira also slightly unclamp the northern part of the Nyiragongo dike surface.A single dike intrusion probably does not induce a normal stress change large enough to cause an eruption at Nyiragongo.However, repeated dike intrusions in Nyamulagira's SSE flank may contribute to the failure of the shallow plumbing system of Nyiragongo.

Influence of the Crustal Extension
The EARS extension probably plays a role in the preferred diking orientations and on the amount of intruded magma volumes.The tectonic rift extension, σt, in this part of the rift is oriented ~N110E [11,12].Resolving this extension on the potential dike directions for host rocks under lithostatic stress,   , we get smaller amplitudes of host rock stresses on the N45E than on the N155E potential dike directions, with values of:   45 =   − sin(110°− 45°)   =   − 0.9   (7) and: respectively.
Next, we assume that failure of the reservoir corresponds to the same magma pressure (Pm) whatever the intrusion orientation, which might be the case around a reservoir, as each intrusion changes the state of stress in the host medium.We find that the overpressure along the N45E surface ∆ 45 =   −   + 0.9   is larger than the overpressure along the N155E surface ∆ 155 =   −   + 0.7   .If we consider that the length of the intrusion is a linear function of the overpressure, because the length is limited either by magma cooling [69] or by the resistance to fracture at the dike tip [70], dike intrusions are expected to be longer along the N45E direction than along the N155E direction.

Magma Storage and Transport
From the detailed geodetic modeling of the 2010 Nyamulagira eruption and the results of [10], we suggest a model of magma transfer and storage associated with a typical upper flank eruption on the SSE flank of Nyamulagira.The onset of the eruption may be marked by over-pressurized magma intruding upward in a dike from a shallow reservoir located at ~4 km beneath the SW caldera depression.Magma propagation is likely driven by the buoyancy of the gas pocket and gas-rich magma at the dike tip, as evidenced by the intense lava fountaining during the first stage of the eruption and the large rv.The magma overpressure keeps the dike open.Just after the eruption onset, the dike may grow laterally due to the decreased buoyancy corresponding to the increased density of the underlying, gas-poor magma, mainly in the direction of the slope, into a weak and fractured area beneath the SSE flank, as observed, for instance, at Piton de La Fournaise [39] and Etna [71].As the overpressure in the reservoir is relaxed, the magma flow rate decreases; magma can no longer reach the summit caldera and erupts only at the lower flank vent.A few days after the eruption onset, the eruption is restricted to a small portion at the base of the fracture.Finally, after a few more days or weeks, the activity stops.
This study provides insights for the shallow magma plumbing system active during 2006-2012 (Figure 7).In 2006, magma intruded from a depth into a shallow reservoir located at an ~1-km depth beneath the SW depression of the summit caldera [20] as a sub-vertical dike beneath the volcano and probably grew laterally in the SSE flank of the volcano, where an eruptive fissure opened.This process of vertical ascent feeding lateral dikes is common for small volcanic edifices (e.g., [39,72,73]).We also suggest that the 2010 sources unclamped the deep part of the potential dike surface oriented ~N45E, south of the caldera, in which magma could have intruded from a deep reservoir.Magma could have then migrated laterally following this ~N45E SC orientation, which may represent a plane of weakness due to pre-existing faults buried under the lava field [13] toward the 2011-2012 eruption site location.Unusual seismic tremor data [10] also suggest magma movement at a depth along a similar ~N45E orientation.

Conclusions
The best-fit model describing the 2010 co-eruptive deformation, which is also the most likely according to the AIC values (Table 1), indicates two subvertical dike intrusions and deflation of a magma reservoir.The northern dike is connected to an intra-caldera ~N45E-trending eruptive fissure, extending to an ~2.5-km maximum depth.The southern dike is connected to an ~N175E-trending flank eruptive fissure, extending to the depth of the inferred reservoir at ~3.5 km.The magma reservoir is located ~3.5 km beneath the SW depression in Nyamulagira caldera.The reservoir location is coincident with the reservoir active during the 1938-1940 [10] and 2006 [20] eruptions.The ratio of total emitted volume to the decrease in reservoir volume decrease is determined to be ~9.3,suggesting that the magma in the reservoir is gas-rich, allowing gas exsolution to restore the pressure loss induced by the eruption.This value is consistent with the shallow and small magma reservoir modeled in this study.
From a normal stress change analysis, we infer that the 2010 intrusions beneath the SSE flank of Nyamulagira encouraged diking at a depth along the ~N155E and ~N45E SC surfaces.The 2010 intrusions might have encouraged the ascent of magma from a deeper reservoir, which then fed an intrusion in 2011-2012.Due to regional rift extension, more extensive intrusions would be favored along the ~N45E orientation, providing a possible explanation for the ensuing distal 2011-2012 eruption.Finally, our normal stress change results also suggest that repeated dike intrusions beneath Nyamulagira's SSE flank slightly encourage intrusions beneath Nyiragongo.

Figure 3 .
Figure 3. Geometry of the preferred best-fit model including two dikes and a deflating reservoir: (a) Map view.The 2010 eruptive fissures are shown as bold green lines.(b) North-vertical cross-section.(c) East-vertical cross-section.Volume changes of the dike associated with the caldera's fissure, the southern flank's fissure and the spherical reservoir are 0.74, 4.2 and −5.1 × 10 6 m 3 , respectively.

Figure 4 .
Figure 4. Close-up of the distribution of recent eruptive fissures at Nyamulagira.The potential dike surface orientations investigated in the normal stress analysis (Section 5.4) and shown in Figure 6 and Figure S3 are represented with brown dashed lines.SC stands for south of the caldera.

Figure 5 .
Figure 5. Sequence of seismic, volcanic and deformation events related to the 2010 Nyamulagira eruption.a-d correspond to the time spanned by the corresponding interferograms shown in Figure2.See[10] for a detailed review of the four eruptive stages.SP and LP stand for short and long period, respectively.

Figure 6 .
Figure 6.Perspective plots from two orthogonal viewing angles showing the changes in the normal stress on the potential dike surfaces (Figure 4) caused by the preferred model of the 2010 best-fitting deformation sources (black mesh).The color scale corresponds to normal stress change in MPa, with positive values (clamping) clipped to dark red and negative values (unclamping) from blue to red.SC stands for south of the caldera.Both the N155 and N45E SC dike surfaces are unclamped over large and deep areas beneath the volcano by the 2010 eruption sources, while only a small area of the N45E (caldera intersecting) surface is unclamped.

Figure 7 .
Figure 7. Schematic shallow magma plumbing system for the period 2006-2012.Arrows represent possible magma migration paths.Red and blue sources experience inflation and deflation, respectively.Orange stars denote the locations of eruptive fissures and vents.SC stands for south of the caldera.The two dikes and magma reservoir inferred from the modeling of the 2010 eruption data are represented with black meshes.

Table 1 .
Comparison of the best-fit model obtained for each combination of sources inverted.Note that several inversions were run for each possible source combination; only the best-fit result obtained for each combination is shown.* denotes the preferred model among all of the best-fit models.Percentage of explained deformation and RMS error are defined in the Supplementary Materials.

Table 3 .
Thermodynamics parameters and major and minor element composition considered for the 2010 lavas of Nyamulagira used to compute the weight percent of the exsolved gases.