The rheology of crystallizing basaltic lavas from Nyiragongo and Nyamuragira volcanoes, D.R.C.

Nyiragongo (Democratic Republic of Congo) erupts fluid, fast-moving, foidite lavas. Nearby Nyamuragira frequently erupts less mobile tephrite lavas. Nyamuragira flows rarely threaten urbanized areas, but Nyiragongo flows threaten Goma (pop. ~900,000) and surrounding villages, resulting in fatalities in 1977 and 2002. We report new laboratory measurements of viscosity evolution during cooling and crystallization from both volcanoes by concentric cylinder viscometry. Melt viscosity is ~33 Pa s at the liquidus (~1220°C) for Nyiragongo lavas, similar to Hawaiian basalts. Lavas remain fluid over ~75°C of undercooling (ϕc 10 ms-1 even at subliquidus temperatures on low slopes, posing great risk to Goma.


Introduction
Lava flows are typically regarded as relatively modest volcanic hazards when compared to other faster and farther-reaching hazards such as pyroclastic density currents, lahars, or debris flows [National Academies of Sciences, Engineering, and Medicine 2017]. Yet, hot, fluid, mafic lavas can advance rapidly and are particularly hazardous when human settlements begin to encroach on the volcanic flanks, for example at Mt Etna, Italy [Guest and Murray 1979], Fogo [Richter et al. 2016], or Kīlauea, USA [Turner et al. 2017]. Nyiragongo volcano, in the Democratic Republic of Congo (DRC), is another example of where lava flows have resulted in considerable societal impacts and also the loss of life [e.g. Allard et al. 2002;Baxter et al. 2002]. Along with its neighbor Nyamuragira volcano, these are the only two active volcanoes in the Virunga Volcanic Province (VVP), part of the East African Rift System (EARS; Figure 1). With the city of Goma located on the southern flank of Nyiragongo, this volcano poses a high risk and was named a "decade volcano" by the United Nations and International Association of Volcanology and Chemistry of Earth's Interior [Newhall 1996] to push research efforts toward understanding volcanoes that Corresponding author: aaron.a.morrison@gmail.com have a high likelihood of impacting human society.
This study contributes to this goal by experimentally investigating the rheology of lavas from these two volcanoes. Rheology is the study of flow, or the way material respond to stress, and links the physical processes occurring during flow emplacement to final post-emplacement morphology. The main factors controlling lava rheology are temperature, composition (± volatiles), and crystal and/or bubble fractions [e.g. Bagdassarov and Pinkerton 2004;Diniega et al. 2013;Mader et al. 2013;McBirney and Murase 1984]. To quantify changes to rheology we look at how parameters such as viscosity and apparent yield strength change as a function of cooling and crystallization. In this study, we conducted super-and subliquidus viscosity experiments to quantify these changes. Such information is essential for the physical modeling of lava flow emplacement and the associated hazards [e.g. Harris and Rowland 2001;Wantim et al. 2013]. Using electron probe microanalysis (EPMA) and scanning electron microscopy (SEM) techniques, subliquidus experiments conducted at different temperatures can be compared to observe changes in crystal content and chemistry of the interstitial melt. This study attempts to quantify experimentally the rheological behavior of lavas from Nyiragongo and Nyamuragira, two volca-2 Background

Previous rheological studies
Many studies have investigated lava rheology in terms of two-phase (melt + crystals or melt + bubbles) and three-phase (melt + crystals + bubbles) materials [Mader et al. 2013;Truby et al. 2015]. Silicate melts behave as Newtonian fluids until they either cool below their liquidus and crystallize a critical volume fraction or are deformed faster than the timescale of relaxation of the melt [Dingwell and Webb 1989]. Concentric cylinder viscometry can be used to measure the viscosity of high temperature (superliquidus) materials [e.g. Bockris et al. 1955;Mysen et al. 1985;Scarfe and Cronin 1986;Urbain et al. 1982]. At much lower temperatures and much higher viscosities, parallel-plate viscometry can be used to measure viscosity just above the glass transition [e.g. Avard and Whittington 2012;Bouhifd et al. 2004;Cordonnier et al. 2009;Villeneuve et al. 2008]. Rheological constitutive equations describe the relationship between applied stress and resulting strain rate. The simplest form of this relationship applicable to lavas is that of a Newtonian fluid: where η is the viscosity andγ is the strain rate. Silicate melts behave as Newtonian fluids except at very high strain rates [Webb and Dingwell 1990]. Measurements can also be made at subliquidus temperatures with the concentric cylinder viscometer, allowing for the quantification of the effect of crystals on viscosity [e.g. Chevrel et al. 2015;Ishibashi and Sato 2007;Ishibashi and Sato 2010;Pinkerton and Norton 1995;Ryerson et al. 1988;Sehlke and Whittington 2015;Sehlke et al. 2014;Soldati et al. 2016;Vona et al. 2011]. As crystals nucleate and grow, a crystal framework may develop as they begin interacting with each other. This framework may begin to resist flow, leading to the appearance of a yield strength (σ y ), the minimum required stress before flow can initiate [e.g. Barnes 1999;Nguyen and Boger 1992]. Once crystals and/or gas bubbles are included, experiments have shown that a linear relationship typically no longer holds [Giordano and Dingwell 2003;Giordano et al. 2008;Sehlke and Whittington 2015;Sehlke et al. 2014]. A more general constitutive relationship is the Herschel-Bulkley equation: where σ y represents the initial yield strength that must be overcome before strain occurs, K is the consistency of the material (equivalent to the viscosity for a strain rate of 1s −1 ), and n is the flow index or power-law exponent, which describes the deviation from linearly proportional (Newtonian) behavior. A flow index greater than one implies shear-thickening behavior, and a flow index less than one implies shear-thinning behavior. Bingham materials are characterized by n = 1 but with a detectable yield strength. Material with no yield strength that still exhibit shear thinning/thickening behavior (n 1) is characterized with a power-law. Giordano et al. [2007] made rheological measurements of Nyiragongo lavas, but they did not quantify how rheology changes with crystal content or strain rate. In the current study, viscosity was measured during cooling experiments until crystallization caused a deviation from Newtonian behavior of the liquid trend. These experiments provide disequilibrium information while the isothermal experiments that we conduct provide complimentary information under equilibrium conditions (i.e. constant temperature and crystal fraction).

Regional Setting
East African Rift volcanism spans the entire range of magmatic compositions, from felsic to ultramafic [Andersen et al. 2012; Barette et al. 2017] The rift system has been suggested to divide into two branches around a microplate centered on the Tanzanian Craton [Calais et al. 2006;Saria et al. 2014]. Within the western branch of this division lies the VVP [Smets et al. 2016], located at the border between the DRC, Uganda, and Rwanda. This province contains eight volcanoes, only two of which are active [Demant et al. 1994]: Nyamuragira (commonly: Nyamulagira) and Nyiragongo. These are the most active volcanoes in Africa [Coppola and Cigolini 2013;Smets et al. 2015;Wright et al. 2015], with Nyiragongo erupting some of the lowest viscosity lavas on Earth [Barette et al. 2017;Chakrabarti et al. 2009a;Chakrabarti et al. 2009b;Giordano et al. 2007]. While these two volcanoes are separated by only 13 km, they produce very different lavas suggesting the volcanoes sample different sources and/or are fed by magma chambers with different degrees of differentiation [Andersen et al. 2012;Barette et al. 2017;Chakrabarti et al. 2009a;Smets et al. 2016].

Nyiragongo volcano
Nyiragongo has a stratovolcanic morphology (Figure 1C) rising approximately 3470 m above sea level, 2000 m above the rift valley plain, and is located 18 km north of Lake Kivu in eastern DRC ]. The volcanic system itself is composed of the main cone of Nyiragongo as well as parasitic cones to the north (Baruta) and south (Shaheru) [Sahama 1978]. This volcano is well known for its persistent lava lake activity [Platz et al. 2004;Smets et al. 2017;Tazieff 1949]. The main crater is roughly 1.3 km in diameter with a lava lake approximately 250 m across making it the largest persistent lava lake in the world [Sawyer et al. 2008;Smets et al. 2017]. Nyiragongo lavas are  Smets et al. [2016], based on previous field observations and mapping from Sahama [1954], Thonnard and Denaeyer [1965], Denaeyer and Schellinck [1965], Smets et al. [2010], and Smets et al. [2015]. Bold black lines are the main rift segments of the western branch of the East African Rift, based on the mapping of Smets et al. [2016]. Thin white lines represent eruptive fissures Smets et al. 2016 [Chakrabarti et al. 2009a;Chakrabarti et al. 2009b] which form thin flows indicative of low viscosity ( Figure 2). Giordano et al. [2007] report viscosity measurements for 2002 Nyiragongo lavas (~60 Pa s at~1100 • C), which have the lowest viscosity of any terrestrial silicate lava, with only carbonatite lavas (~0.3-5 Pa s) being more fluid [Dawson et al. 1990].
While some studies disagree about the timing of early activity of this volcano, there is consensus that persistent lava lake activity occurred at least from 1928 until the 1977 eruption that drained the caldera [Demant et al. 1994;Simkin and Siebert 1994;Tazieff 1984;Von Götzen 1895]. There have been only two major flank eruptions in the past 50 years (1977 and 2002) that resulted in fatalities in the nearby city of Goma and its vicinity. A roughly north-south oriented fracture system propagated southward from the main summit crater allowing lavas to infiltrate the city of Goma [Acocella and Neri 2009;Carn et al. 2003]. In 1977, the city had a population of around 70,000 [Barigora 2008]. Lava lake activity returned for a few months in 1982[Durieux 2002aKomorowski et al. 2002;Krafft and Krafft 1983;Tazieff 1984;Tedesco 2003].
Subsequently, the city expanded westward and northward toward the volcano, growing in population to over 400,000 people in 2002 [Barigora 2008], largely due to the influx of refugees from the recurrent armed conflicts in this region. This set the stage for the eruption in January of 2002 to be one of the most devastating eruptions in recent history. A very thorough chronology of the 2002 eruption is compiled in Komorowski et al. [2002]. During this eruption, the fracture network from 1977 reactivated and extended southward entering the city of Goma. The extremely low viscosity lavas rapidly infiltrated the town via these fissures [Carn et al. 2003;Komorowski et al. 2002;Sawyer et al. 2008], resulting in 50 deaths directly caused by the lava itself and displacing almost the entire population [Allard et al. 2002;Baxter et al. 2002;Komorowski et al. 2002]. Lavas that erupted at the foot of the main volcanic edifice were relatively gasrich, driving fountaining processes, in contrast to the relatively degassed lavas that escaped from the fissures located between the Shaheru crater and the Nyiragongo main cone [Giordano et al. 2007].

Nyamuragira volcano
Nyamuragira is a shield volcano ( Figure 1C), rising approximately 3058 m above sea level, located 25 km north of Lake Kivu in eastern DRC and only 13 km northwest of Nyiragongo [Smets et al. 2014b]. This volcano has been erupting approximately once every 1-4 years, with 38 reported eruptions since the 1900s . Generally, lavas from Nyamuragira are alkali basalts, hawaiites, basanites and tephrites ( Figure 2) with a range of silica content from 43-56 wt.% [Aoki et al. 1985;Barette et al. 2017].
Over the last 25 years, Nyamuragira has been producing flank eruptions every few years with activity in : 1991-1993, 1996, 1998, 2000, 2001, 2002[Coppola and Cigolini 2013Smets et al. 2015;Wadge and Burt 2011, and references therein]. In 2014, intermittent summit activity developed in the northeastern part of the main caldera, producing effusive or lava lake activity [Smets et al. 2014a]. Most flank eruptions at Nyamuragira have a very abrupt and vigorous beginning that can be very well constrained temporally. However, activity at the end of an eruption diminishes gradually making an end date much more difficult to establish in the absence of daily field observation . These regular eruptions are probably fed by one or more shallow, pressurized magma chambers, which source magmas from depth [Burt et al. 1994;Smets et al. 2015]. Smets et al. [2015] characterized four different categories of eruptions. The first are summit and upper flank eruptions that reflect changes in the fractured edifice. The second are classical flank eruptions that reflect the fissure system crossing the gravitational stress   (NYI-2002 andNYA-1948) in relation to some of the other mafic compositions that have experimental measurements of rheology (Hawaii: Sehlke et al. [2014]; Iceland: Kolzenburg et al. [2017]; Pacaya: Soldati et al. [2016]; Alkaline Olivine Basalt: Ishibashi and Sato [2007]; Picrite: Ryerson et al. [1988]; Fuji: Ishibashi and Sato [2010]; Stromboli and Etna: Vona et al. [2011]). The range of compositions for both Nyiragongo lavas and Nyamuragira lavas, from Barette et al. [2017], are shown as blue and orange fields respectively. field. These first two categories are eruptions that likely source shallow magma chambers at 3-7 km depth. The third are long lived flank eruptions (>65 days) which may source deeper magma chambers [e.g. Wauthier et al. 2013] but have other varying interpretations. The final category are atypical eruptions that do not fall into the other three categories and seem to be unconnected to the main central plumbing system of Nyamuragira. near the flow front. The Nyamuragira samples were taken approximately 20 km from the flow front near the edge of Lake Kivu, west of Goma. Petrographic thin sections were made for each composition. Images of petrographic thin sections were used to estimate textural characteristics (see subsection 4.8). The NYI-2002 sample was~45 % phenocrysts which were dominantly plagioclase. The matrix was~45 % and was composed of plagioclase, nepheline, and pyroxene microlites. Vesicularity was~10 %. The NYI-1977 sample was~10 % plagioclase and nepheline phenocrysts. The matrix was~80 % of the section composed again of plagioclase and nepheline microlites. Vesicularity was 10 %. The NYA-1948 sample was~25 % phenocrysts of mostly plagioclase. The matrix was~60 % of dominantly plagioclase and pyroxene microlites. Vesicularity was~15 %. NYI-1977 bulk composition is within error of NYI-2002 apart from Fe and K (Table 1). Phenocrysts were subhedral to euhedral and no glass was observed in any of the natural samples. 4 Experimental methods

Glass synthesis
Samples were individually crushed in a steel shatter box for two minutes. The powder was then collected in a previously iron saturated Pt 90 -Rh 10 crucible and placed in a muffle furnace where it was heated to 1600°C. To ensure homogenization, the melt was swirled regularly during the heating process. After two hours the melt was quenched on a copper plate and allowed to cool before remelting and requenching a second time to aid in homogenization.

Differential scanning calorimetry
Liquidus temperatures were determined using a Netzsch DSC 404 Pegasus F1 differential scanning calorimeter (DSC) with Pt-furnace. A previously ironsaturated Pt-Rh alloy was used as the sample and reference container material. Each measurement consisted of a 15-minute isothermal equilibration segment (typically 50°C), a 20 K/ min heating segment up to 1500°C , a 5-minute isothermal equilibration segment, and an uncontrolled cooling segment where the furnace was turned off. Liquidus temperatures were determined as the point of final melting from the heat flow curves (i.e. the end of the final endothermic peak). These values are reported in Table 1.

Parallel plate viscometry
Cylindrical cores (approximately 6 mm diameter by 8 mm height) were drilled from the remelted, crystalfree, glass samples. A Theta Industries Rheotronic III parallel-plate viscometer was used to measure viscosity between 600 and 800°C under uniaxial compression applied by a 1500 g load. This instrument is capable of viscosity measurements between~10 8 and 10 13 Pa s. A transducer is used to measure sample length to ± 0.1 µm and viscosity (η) is calculated using the following relationship [Dingwell 1995], assuming constant sample volume and perfect slip: In this equation, m is the mass of the applied load, g is gravitational acceleration, h is the sample height, V is the sample volume, and (δh/δt) is the rate of change in height. The viscometer was calibrated using National Institute of Standards and Technology (NIST) standard glasses and verified to be accurate to within 0.06 log units, with a precision based on repeat measurements of better than 0.04 log units [Whittington et al. 2009]. Viscosity data for the crystal-free supercooled liquid of the 2002 Nyiragongo flow were previously collected in Sehlke and Whittington [2016].

Concentric cylinder viscometry
A Theta Industries Rheotronic II Rotating Viscometer was used with a Brookfield DV3TRV Rheometer head for superliquidus experiments. This head has a spring torque value calibrated to 7.187 × 10 -4 N m which, in our setup, is capable of measuring viscosities between 0.1 and 1000 Pa s. Instrument calibration was checked against Brookfield standard oils ranging from 10 to 100 Pa s. Remelted glass shards were placed in a Pt 90 -Rh 10 crucible 65 mm tall and with an internal radius of 16.075 mm. The crucible is then held by three alumina rods and lowered into a box furnace capable of operating to 1600°C. A Pt 90 -Rh 10 sleeve with a hemispherical tip is attached to an alumina spindle by a through-going Pt wire, which is then lowered into the melt. The sleeve has a length of 65 mm and an external radius of 3.75 mm. A particular rotation rate is chosen, and the torque required to achieve that rate is measured. Shear stress (σ ) is then defined by: where τ is the torque on the spindle, R b is the radius of the spindle, and L is the immersion depth of the spindle. Immersion depth is controlled by a micrometer and actual depth is measured directly by residual melt left on the Pt 90 -Rh 10 sleeve after each experiment. With this setup, the assembly has a wide gap geometry such that the ratio of crucible radius (R c ) to spindle radius (R b ) is approximately 4.02. This allows crystallizing phases with large aspect ratios, such as plagioclase, to align with the flow direction. The following equation is used to calculate shear strain rate (γ): where ω is the angular velocity of the spindle and n is the flow index calculated by finding the slope of a linear regression of ln(τ) and ln(ω) [Stein and Spera 1998]. Superliquidus viscosity data for the 2002 Nyiragongo flow has been previously collected in Sehlke and Whittington [2016]. Parallel plate and concentric cylinder data were then fit using a Vogel-Fulcher-Tammann (VFT) equation [Vogel 1921].
Subliquidus rheological experiments were performed at various crystal fractions using a Brookfield HBDV-III Ultra measuring head with a larger spring constant of 5.7496 × 10 −3 N m. NIST standard glasses were used to check the accuracy of the viscometer, which was determined to be better than ± 0.06 log units [Getson and Whittington 2007]. Samples were heated to 1450°C at a rate of 5°C min −1 to ensure complete melting of all crystal phases. Viscosity was verified at this temperature for each run to ensure consistency across experiments. The spindle was then inserted and allowed to rotate at a constant angular velocity while the sample cooled at 10°C min −1 to a target temperature below the liquidus. Once the target temperature was reached it was allowed to equilibrate overnight for 10-12 hours, until the torque stabilized for at least four hours. After the equilibration period, the strain rate was varied by changing the angular velocity at regular intervals to measure over the widest range of torque values as possible, between 10 to 80 % of total spring capacity. Each segment at a given angular velocity lasted for 45 to 60 minutes to allow for the hysteresis effects to be overcome. Data from the last 10 to 20 minutes of each segment where torque readings were stable were used in calculations. After all programmed segments were completed, the spindle was removed and the crucible quenched immediately in a room-temperature water bath for 5-10 minutes to prevent further crystallization during cooling. After reaching room temperature, glass cores were drilled from the crucible and sampled at varying depths from the top to approximately 3 cm depth. The interstitial glass composition of the lowest subliquidus experiment for each composition was then synthesized from reagent grade oxide powders and measured following the same procedure as the initial liquid measurements.

Density determination
Densities were measured after each experiment using the Archimedean method. Core sample weight was measured both in air (w air ) and while immersed in anhydrous ethanol w eth . Density of the sample can then be calculated as follows: where ρ eth is the temperature corrected density of the ethanol. Five measurements were made on each sample as well as a standard quartz reference with sample precision of ± 8.4 kg m −3 (2σ ), or better than 0.4%. Data are reported in Table 2.

Chemical analysis
The remelted glass composition was analyzed by EPMA with results summarized in Table 1. Nyiragongo lavas plot as foidites and Nyamuragira lavas plots as tephrites ( Figure 2). Each composition was checked to verify that the glass was crystal free and composition was homogenous from at least two different chips of glass. Chemical compositions of interstitial glass and individual crystals from experiments were analyzed using a JEOL JXA-8200 electron microprobe, equipped with five wavelength-dispersive spectrometers and a JEOL (e2v/Gresham) silicon-drift energydispersive spectrometer, located at Washington University in St. Louis. The apparatus was operating with a 15 kV accelerating potential and 25 nA probe current for all spot analyses. A 5 µm diameter electron beam was used for crystals to avoid volatilization of sodium and a 20 µm diameter electron beam was used for glasses. Silicate and oxide standards were used for calibration and data were corrected using CITZAF following Armstrong [1995]. Approximately ten individual analyses were obtained for glasses and three for each crystal phase present.

Fe-redox
Long-duration subliquidus experiments in air could potentially result in a change in sample oxidation state, which may have an effect on sample viscosity [e.g. Dingwell and Virgo 1987]. We determined the redox state of Fe in the natural rock samples, the remelted bulk glasses, and the quenched products of each subliquidus experiment. Bulk iron oxidation state was measured using colorimetric methods described by Wilson [1960] and Schuessler et al. [2008]. Ultraviolet/Visible (UV/Vis) photo-spectrometry was carried out with a Nanodrop1000 UV/Vis spectrometer. Accuracy was checked against a blank and USGS standards (BIR-1a basalt and W-2a Diabase) to an uncertainty of 0.20 wt.% (2σ ). A full description of our procedure can be found in Sehlke et al. [2014].

Crystal volume fraction and aspect ratio
For each subliquidus viscosity experiment, quenched samples were drilled and 2-4 pieces recovered for crystal volume fraction calculations. As the samples shattered when drilled, the exact location of the recovered sample from within the crucible cannot be precisely determined, but, qualitative location (i.e. top, middle, bottom) within the core is known. Several backscattered electron (BSE) images were acquired on each piece at magnifications of 40×, 180×, and 2000×. These images were processed using Adobe Photoshop® to isolate single phases present in the sample. The number of pixels were counted for each phase using the histogram function and the average of results using three different thresholds (or shading tolerances) is reported. This average is to account for the inexact crystal boundaries formed by the square pixels. This inexactitude may result in melt and crystal fractions that do not sum to 1.0 but are still within a 2σ uncertainty. Processing images at different magnifications for the same sample results in a volume fraction for each phase that has a standard deviation of 2-4 vol.%. Approximately 6-10 images were used for each sample piece. This method was also used for petrographic thin sections to calculate vesicle, phenocryst, and matrix fractions. The freeware program ImageJ (http://imagej.nih.gov/ij/) was used to determine crystal aspect ratios by utilizing the automatic ellipse fitting procedure to the individual grain boundaries for each present phase. Each ellipse was subsequently checked for accuracy manually.
If the program were unable to adequately represent the aspect ratio with an ellipse, aspect ratio was calculated by hand. The mean aspect ratio (R) for a given experiment is calculated using the following equation: where R i is the apparent mean aspect ratio of each phase and φ c is the crystal area fraction.

Calorimetry
The liquidus temperature for Nyiragongo was determined from calorimetry, during heating of lava samples collected in the field, to be~1220°C (Table 1). This is supported by SEM images (Figure 3) of the 1221°C experiment where only a few individual crystals were observed across six different SEM images. Based on calorimetry, the liquidus temperature for Nyamuragira is~1260°C. This is consistent with our rheology experiments, where we observe a minuscule crystal fraction (φ c = 0.001 ± 0.004) determined from SEM images (Figure 4).

Liquid viscosity
Compositional data of the quenched glasses analyzed by EPMA are summarized in ). These data are presented in Table 3, fitted with a VFT equation [Vogel 1921], and plotted on an Arrhenian diagram in Figure 5. The Nyamuragira liquid is 1-2 orders of magnitude more viscous than the Nyiragongo liquid: −0.20 log Pa s at 1550°C and 10.50 log Pa s at 700°C for NYI-2002 compared to 0.23 log Pa s at 1550°C and 11.78 log Pa s at 700°C for NYA-1948.   (NYI-2002). The black scale bar in each image is 500 µm. Crystal fraction uncertainty is roughly ± 0.02 (see Table 2). White crystals are spinels (sp) and light gray crystals are clinopyroxene (cpx). Figure 4: Backscattered electron images of each subliquidus experiment for Nyamuragira (NYA-1948). The black scale bar in each image is 500 µm. Crystal fraction uncertainty is roughly ± 0.02 (see Table 2). White crystals are spinels (sp) and dark gray crystals are plagioclase (pl).  Figure 5: Arrhenian diagram plotting viscosities for the bulk remelted material. Circles and diamonds represent measured Nyiragongo and Nyamuragira data respectively. Solid curves are data from Table 3 fit using models from Vogel [1921]. Dotted lines are predictions of the model from Giordano et al. [2008] and dashed lines are predictions of the model from Sehlke and Whittington [2016].

Subliquidus experiments
Eight subliquidus experiments were conducted for NYI-2002 between 1221-1145°C. Eight subliquidus experiments were also conducted for the NYA-1948 between 1255-1151°C. BSE images for each subliquidus experiment are presented in Figure 3 and 4 for both Nyiragongo and Nyamuragira. Table 2 compiles crystal fraction, aspect ratio, and mineralogical data for each image displayed. The first NYI-2002 experiment was done at the approximate liquidus temperature of 1221°C . No significant crystallization occurs in the NYI samples until~50-60°C undercooled. Oxidizing conditions suppress olivine growth and facilitate crystallization of Fe-Mg spinels (white). Augite (light gray) begins crystallizing at 1150°C but not abundantly within the experimental temperature range. The highest crystal fraction achieved was 0.04 at 1145°C. The first NYA-1948 experiment was done at the higher near-liquidus temperature of 1255°C. Crystallization occurs in the Nyamuragira samples in a much smaller range of un-dercooling~20-30°C. Fe-Mg spinels again are the first phase to crystallize with plagioclase (dark gray) coming in at 70-80°C undercooled. The highest crystal fraction achieved was 0.11 at 1151°C.

Melt chemistry and evolution
Crystal fraction increases, melt fraction decreases, and the residual liquid composition becomes more silicic and less iron-rich as temperature decreases ( Figure 6). The iron content decreases more steeply at first as spinel is the only phase forming but decreases more shallowly at lower temperatures as a second phase begins to crystallize. Spinels become slightly more magnesium-rich as temperature decreases. The decrease in Mg seen in the Nyiragongo samples is due to the crystallization of clinopyroxene at lower temperatures. In the CaO curves for each composition, a slight increase followed by a slight decrease is observed. The increase is due to spinel formation which does not incorporate much Ca, effectively enriching the melt. The subsequent decrease is due to the formation of augite in the lower temperature Nyiragongo samples and a calcium-rich plagioclase phase in the Nyamuragira samples. The alkalis also show a slight increase with decreasing temperature as they are incompatible in the crystallizing phases. Compositional trends for Nyiragongo and Nyamuragira are broadly similar over the range of experimental temperatures. All oxidation state information is reported in Tables 1 Figure 6: Melt chemistry evolution for Nyiragongo (blue circles) and Nyamuragira (orange diamonds) compositions. Data points represent the average of 6-8 EPMA analyses of interstitial glass from a given subliquidus temperature. Onset of pyroxene crystallization for the NYI sample and the onset of plagioclase crystallization for the NYA sample are indicated by the arrows. and 4. The Nyiragongo bulk compositions (NYI-2002 andNYI-1977) show a marked dichotomy in oxidation state. Both samples have similar total FeO wt.% as measured by EPMA and UV/Vis spectroscopy varying by only~0.50 wt.% (Table 1). However, the sample from the 1977 eruption is significantly more oxidized compared to the 2002 eruption (Table 1). Nyamuragira has a similarly high iron content at 11.62 ± 0.10 wt.% and shares a similar oxidation state as the 2002 Nyiragongo eruption. The quenched subliquidus samples are all more similar in oxidation state to the bulk NYI-1977 material (Table 1) than the more reduced bulk NYI-2002 material. Table 4 compiles the results of these calculations for each subliquidus experiment.
The flow index (n) was determined from the slope of a linear regression of ln(τ) versus ln(ω) for each subliquidus experiment. Linear regression parameters and 2σ uncertainties are reported in Table 2. Both compositions behave as Newtonian fluids (n = 0.98 ± 0.01) above the liquidus. Nyiragongo continues to exhibit near-Newtonian behavior (flow indices between 0.94-0.96) at nearly 50°C below the liquidus temperature due to the low amount of crystallization over that temperature range. Flow indices decrease at lower temperatures and higher crystal fractions. At the lowest temperature of 1145°C, corresponding to~60°C of undercooling, the crystal fraction is 0.04 ± 0.02 and the flow index is 0.80 ± 0.06. Bulk viscosity for Nyiragongo increases from 35 Pa s (at~6 s −1 ) at the liquidus to 644 Pa s (at~3 s −1 ) at the highest observed crystal fraction.
Crystal fraction begins increasing much sooner upon undercooling in the Nyamuragira samples. At~25°C below the liquidus the first significant crystallization occurs where Fe-Mg spinels account for a crystal fraction of 0.02 ± 0.02. Non-Newtonian behavior (n ≤ 0.95) is observed almost immediately (1255°C, n = 0.89 ± 0.02) upon crossing the liquidus. Flow indices remain constant to within uncertainty of each other over~100°C until the lowest temperature experiment at 1151°C where n drops to 0.48 ± 0.12. This decrease is attributed to the much higher crystal fraction (0.11 ± 0.02; about twice that of the highest Nyiragongo sample) with this lava before the limitations of the instrument were reached. Bulk viscosity for Nyamuragira increases from 38 Pa s (at~4 s −1 ) at the liquidus to 559 Pa s (at~4 s −1 ) at the highest crystal fraction.
Several attempts at temperatures below~1145°C and~1151°C for Nyiragongo and Nyamuragira respectively, produced unreliable data. At these temperatures, a surface crust (~0.5-1 mm thick) formed a more rigid layer that did not flow in the same manner as the melt and thus, the measurement ceased to accurately reflect the sample viscosity. The presence of surface lobate structures show that this crust was still able to flow but it was no longer homogenous with the rest of the melt. This is likely the result of oxidation at the interface between the molten sample and the air.

Yield strength
Once a crystal network begins to form, a minimum stress (yield strength) may be required to overcome the resistance of this network to flow. This yield strength (σ y ) was estimated as the nonzero, y-intercept of the flow curve (Figure 7) on the stress axis [Barnes 1999;Bingham and Green 1919;Bingham 1922;Kerr and Lister 1991]. Due to the finite sensitivity of the instrument and infinite time required to measure at a strain rate of 0 s −1 , data from the flow curves are linearly extrapolated back to the stress axis to estimate a yield strength. Because this linear fit represents a Bingham rheological behavior, these estimates reflect the maximum possible value of apparent yield strength for each experiment.
For the Nyiragongo lavas, Bingham yield strength estimates are either within error of zero or merely a few tens of Pa (Figure 8). Similarly, the results for Nyamuragira lava are mostly within error of zero. The highest crystal fraction experiment (1151°C; φ c = 0.11 ± 0.02) resulted in a yield strength estimate of 490 ± 101 Pa. All yield strength estimates are reported in Table 2 and Figure 8.

Estimation of Nyiragongo liquidus temperature
We determined the liquidus temperature to bẽ 1220°C by calorimetry and by visual confirmation of crystals in quenched experimental products (Figure 3). Giordano et al. [2007] made rheological measurements on similar lava samples erupted from the 2002 fissure network. They conducted dynamic cooling experiments (1, 2, 3, 5 K min −1 ) in which they estimated the liquidus temperature by the departure from the pure liquid viscosity trend (~1100°C). This method always underestimates the true liquidus temperature due to undercooling, which becomes more pronounced at higher cooling rates ( Figure 9). Furthermore, even in our experiments, the apparent viscosity closely follows the trend of the bulk liquid until~1165°C, after more than 50°C of cooling below the liquidus, because only at this temperature does the crystal fraction become large enough to detect from rheology measurements ( Figure 10). We have plotted our data alongside that of Giordano et al. [2007] to illustrate the difference between cooling and isothermal experiments (Figure 9). The cooling paths of Giordano et al. [2007] show the viscosity trends as samples cool at different (constant) rates. Our data, collected from isothermal segments, represent the equilibrium condition which is effectively the 0 K min −1 cooling experiment. Our data thus create the end-member of those of Giordano et al. [2007], which is useful for putting realistic bounds on the possible achievable viscosity range.
Ultimately, neither data set perfectly replicates nature since lavas do cool as they move (unlike isother-    Bingham (open) and Herschel-Bulkley (closed) fits were used and values were also calculated following Dragoni [1989] (colored). Error bars for these calculated values are smaller than the symbol size.
mal measurements) but not at a constant rate (unlike cooling experiments). However, both data sets are complimentary and provide a better understanding of what the full range of possible viscosity paths could be. We suggest that both types of experiments should be run in future rheological studies.

Crystallization
The Rhyolite-MELTS program [Ghiorso and Sack 1995;Gualda et al. 2012] was used to compare the observed crystallization characteristics (crystal fraction and phase compositions) with those modeled by the program. When using the oxygen fugacity inferred from the measured bulk Fe 2+ /Fe 3+ ratios, the model calculates liquidus temperatures that are unrealistic. This suggests that the samples are initially governed by their intrinsic, more reduced oxygen fugacity and become more oxidized during the experiment as diffusion of oxygen into the samples occurs over the experimental timescale. We adjusted the redox conditions of the simulations in order to replicate the liquidus temperatures determined from calorimetry. The resulting calculated log f O 2 values of ∆NNO-0.63 and ∆NNO+0.75 for Nyiragongo and Nyamuragira respectively, also result in the model replicating our experimental crystal fractions and compositions very closely. Figure 11 shows the plots of calculated and observed cooling histories under these internally-buffered f O 2 Figure 9: Nyiragongo viscosity data from Giordano et al. [2007] compared to our data from Figure 10. Blue points are the average apparent viscosities for each subliquidus experiment.
conditions. Rhyolite-MELTS accurately predicts what is observed in the experimental Nyiragongo samples. There is very little crystallization that takes place until at least 50°C of undercooling. Crystal fractions of about 0.04 are achieved prior to the rheological measurement limit, which is clearly seen on these graphs as a sharp increase in slope. This rise in crystallization rate corresponds to a rapid rise in torque, and thus viscosity, which quickly exceeds the limits of the instrument. The break in slope effectively represents the lower temperature limit of experiments and explains the many failed attempts at measuring rheology at lower temperatures. This point also represents where the lava flows can be expected to quickly cease moving as crystal fraction, viscosity, and yield strength quickly increase. As seen in Figure 9, the temperature where such a slope break occurs can be pushed to lower temperatures by increasing the cooling rate of the lava. Our data provide a maximum temperature corresponding to the equilibrium case. Prior to the break in slope, an iron-magnesium spinel is the only phase that forms, both in the lab and in the model. The experimental samples match the calculated curves, terminating at the rapid influx of clinopyroxene where continued measurement was impractical. Augite crystals also are in good agreement compositionally between the observed and modeled values.
The experimental Nyamuragira samples also closely followed the modeled crystallization trend, with an initial gradual increase followed by a rapid increase due to the appearance of additional phases at about 100°C of undercooling. Again, iron-magnesium spinels are the first phase to form and show agreement in composition between observed and calculated values. However, the second phase predicted to appear would be clinopy-A B Figure 10: Viscosity modeled for the bulk remelted material (light line) and the most evolved subliquidus melt composition (dark line). Points plotted are apparent viscosities of each successive subliquidus experiment.
[A] Nyiragongo (blue) shows little crystallization upon cooling as the apparent viscosities follow the bulk material until the~40°C undercooled. Only after significant crystallization causes deviation from the bulk trend can a linear interpolation of interstitial melt viscosity be applied. [B] Nyamuragira (orange) shows a more linear crystallization trend across the range of experimental temperatures. roxene similar in composition to those found in Nyiragongo, but no augite crystals were observed in these experimental subliquidus samples. Instead, plagioclase was the next phase observed while Rhyolite-MELTS predicted the first plagioclase phase to form~15°C after the appearance of augite .
With such a silica-poor, depolymerized composition as Nyiragongo, rapid crystallization might be expected. Taking into account the duration of the experiment, a 10-12-hour equilibration period where torque becomes stable at the target temperature followed by 10-12 hours of isothermal measurement, we assume that an equilibrium crystal fraction is achieved. Since crystallizing phases are either more dense (spinel, augite) or Volcanica 3(1): 1 -28. doi: 10.30909/vol.03.01.0128 less dense (plagioclase) than the liquid, the possibility of crystal sinking or floatation must be considered. This can be checked by mass balance (Table 2): if the bulk assemblage (i.e. composition of the crystals weighted by crystal fraction and composition of interstitial glass weighted by glass fraction) at the end of experiments does not match the starting bulk composition, crystal fractionation may have occurred. A linearized least squares approach was used to determine what fraction of the crystal and interstitial glass compositions are required to reproduce the bulk, starting glass composition. Results of the mass balance based on crystal fractions observed in BSE images indicate a good match with the starting bulk composition, implying that little to no settling occurred on the timescale of the experiment. A Stokes' Law calculation (v = 2∆ρgR 2 /9η) was also conducted with the most favorable viscosity (30 Pas), density contrast (800 kg m −3 ), and crystal size (50 µm), which suggested maximum possible settling distances of approximately a centimeter over the timescale of the experiment. However, we do not see crystal sizes this large or viscosities this low at the same time and thus we calculate settling times of less than a millimeter for crystal sizes and viscosities appropriate for any given experiment. Continuous stirring of the sample also hinders crystal settling, promoting a homogenous distribution of crystals.

Effect of redox state on sample properties
Iron may have different roles in the crystal structure depending on oxidation state. The usual assumption is that ferrous iron (Fe 2+ ) acts as a network modifier while ferric iron (Fe 3+ ) can act as a network former, although both ions can potentially adopt either role (see review by Mysen and Richet [2018]). The different coordination numbers adopted by iron atoms will affect the magma viscosity differently, i.e. higher viscosity with more Fe 3+ and lower viscosity with more Fe 2+ . Therefore, oxidation state should be considered when determining how viscosity changes as the melt chemistry evolves. However, Chevrel et al. [2014] found that the difference between liquid viscosities for alkalic iron-rich basalts measured in air and at low oxygen fugacities (∆QFM-3) were only between 0.01 and 0.02 log units. These compositions are rather similar to those of the Nyiragongo and Nyamuragira lavas, and thus redox state is not expected to affect the melt viscosity significantly in these experiments. It may, however, affect both the temperature of crystallization and the crystallizing assemblage. This is likely the reason why samples in the lab do not perfectly reproduce natural sample mineralogy or crystallization path.

Magma flow laws and yield strength
The calculated yield strengths approach zero, and flow indices approach one, indicating Newtonian behavior A B Figure 11: Crystallization histories for Nyiragongo (circles) and Nyamuragira (diamonds). The colored lines are the data calculated for NYI-2002 (blue), andNYA-1948 (orange) bulk material using the Rhyolite-MELTS program [Gualda et al. 2012]. Curves are constructed assuming an internally buffered f O 2 of −0.63 and +0.75 NNO for Nyiragongo and Nyamuragira respectively. The black lines are the observed crystallization data from the subliquidus experiments.
[A] Data plotted against temperature for modeled and observed values.
[B] Data plotted against the amount of undercooling for a more direct comparison. over~60°C of undercooling for the Nyiragongo experiments (Table 2). Non-Newtonian behavior was only observed for the lowest temperature (1145°C) experiment, which had a flow index of 0.80 ± 0.06 and an apparent yield strength of 163 ± 126 Pa (Figure 8). The results for Nyamuragira are qualitatively similar, however, non-Newtonian behavior is observed at all subliquidus temperatures as the flow index is always lower than for Nyiragongo and less than one. For the lowest temperature experiment (1151°C), flow index (n) drops abruptly to 0.48 ± 0.12 ( Table 2). The large er-ror bars arise from the very limited range of strain rates that could be achieved at the lowest temperature. Similarly, the calculated yield strengths remained within uncertainty of zero until the lowest temperature experiment, where 490 ± 101 Pa (Figure 8) is the first likely real value. Using root mean square deviation (RMSD) as the indicator of fit quality, Bingham and Herschel-Bulkley equations show a better fit than Newtonian or power-law relationships. Considering the lack of crystallization upon cooling and the negligible yield strengths, the Newtonian viscosity should be acceptable for approximately the first 50 • C of undercooling. A Herschel-Bulkley equation is the most rigorous equation to define rheology since it has enough degrees of freedom to account for both shear-thinning behavior and yield strength development while also allowing for the other equations (Newtonian, Bingham, power-law) to be recovered. However, when modeling, it is typically advantageous to reduce the number of parameters for ease of computation. Thus, a Bingham rheology still has a good statistical fit to the data and can be used to approximate viscosity for the range of viscosity over which we were able to measure. A Herschel-Bulkley equation should be used exclusively for extrapolations to higher crystal fractions and lower temperatures where both shear-thinning behavior and yield strength development would be pronounced.
Lavas from both volcanoes remain very fluid for at least 80°C below their liquidus temperatures. Liquid viscosities for Nyiragongo lavas are 0.5-1 order of magnitude lower than those for Nyamuragira (Figure 5). Crystallization occurs at much lower temperatures for Nyiragongo lavas (Figure 11) so apparent viscosities and yield strengths for a given temperature will be drastically different between the compositions (Figure 8).
6.5 Physical effect of crystals on relative viscosity of magma By decreasing temperature in successive experiments, it is possible to determine how increasing crystal fraction affects rheological behavior, as long as the effects of changing residual melt composition are also accounted for [Getson and Whittington 2007]. Relative viscosity (η r ) is defined as: where η s is the apparent viscosity of the suspension and η l is the viscosity of the coexisting, interstitial liquid at that temperature. The liquid viscosity can be determined from direct measurement or by calculation. As discussed in Sehlke et al. [2014] and Sehlke and Whittington [2015], the direct measurement of a synthesized interstitial melt composition is a much more reliable method than calculating viscosities using models such as Hui and Zhang [2007] or Giordano et al. [2008].
The interstitial melt compositions of the lowest temperature, highest crystal content, subliquidus experiments of both lavas were determined from EPMA and then synthesized using dry, oxide powder reagents. Viscosity data for both synthesized, evolved compositions (eNYI and eNYA) are listed in Table 3. In the case of Nyiragongo, the lack of crystallization over the first 50°C of undercooling causes the relative viscosity to track along the bulk viscosity trend. Figure 10 shows the modeled viscosity of the initial bulk liquid and the evolved melt composition with measured apparent viscosity data from each subliquidus experiment. It is inferred that the effect of crystals on the relative viscosity is ultimately negligible over this temperature range.

Viscous heating
Viscous heating may also occur when some of the stress driving fluid flow is dissipated by conversion to heat. For the case of lava (a viscous fluid), the heating rate is calculated as follows: where η is the viscosity and ∂ε/∂t is the strain rate. The viscous heating component has been numerically modeled for a range of lava compositions and helps explain the onset of non-Newtonian behavior observed by some studies [Avard and Whittington 2012;Cordonnier et al. 2012;Costa and Macedonio 2005a;Dragoni et al. 2002;Hale et al. 2007]. This viscous heating effect may play an important role in channel flow dynamics where viscous friction can locally raise the flow temperature along the channel margins resulting in a viscosity decrease and velocity increase [Costa and Macedonio 2005a]. This effect is critical in the context of Nyiragongo where flows already remain fluid for an extended range of undercooling. With such low viscosity flows, shear rates are expected to be very high for these thin flows. The additional heat could allow for these already fast-moving flows to extend further from the vent. The Grätz number (Gz) is a dimensionless number that is useful for describing the temperature distribution across a lava flow (Pinkerton and Wilson, 1994). It is calculated as follows: where E is the effusion rate, d is the flow thickness, κ is the thermal diffusivity, x is the flow length (distance from vent), and w is the flow width. Typically, a Grätz number of~30 represents a flow that has decreased 1 % of the difference between eruption temperature and ambient temperature. Pinkerton and Sparks [1976] looked at cooling limited flows on Etna and found that advancement was arrested when the Grätz number reached~300. Since this number is dependent on the flow length, we can attempt to use the Grätz number to determine whether viscous heating may have extended the range of the flow. An effusion rate of 150 m s −1 was taken from Komorowski et al. [2002] for the 2002 eruption of Nyiragongo and flow thicknesses from our own field observations are between 0.5 m and 2 m. The thermal diffusivity of NYI-2002 sample was previously measured by Hofmeister et al. [2016] and we have taken the value of 3.5 × 10 −7 m 2 s −1 for our calculation. Flow widths were estimated from the maps provided in Chirico et al. [2009], Favalli et al. [2009], and Smets et al. [2010]. For Nyiragongo, flow length estimates range from 2.6-10 km for the 1977 and 2002 eruptions according to the flow mapping of Smets et al. [2010]. In two cases for the 2002 eruption, flow lengths were limited by other factors than cooling or volume. One flow was limited by Lake Kivu and another limited by the filling of the Baruta construct. These flows were not used in Grätz number calculations. The remaining flow lengths are typically the same as or longer than lengths calculated by assuming a Grätz number of 300. This suggests that viscous heating may contribute to the transport and emplacement of lavas from Nyiragongo, however, only at extreme conditions of high speeds (>10 m s −1 ) and small thicknesses (<1 m). It is important to note that the uncertainty on the input variables is very high and the results of this calculation are not fully conclusive. The effect of bubbles on viscosity has also been neglected however and may also play a role in the flow extent.
6.7 Implications for volcanic hazards Favalli et al. [2009] attempted to model future flow paths from different vent locations around the city of Goma to analyze susceptibility to lava flow invasion. However, no attempt was made to incorporate rheological properties. While this may not be a problem for probabilistic models that assume a volume limited flow, rheology does play a critical role in determining lava flow speeds and where cooling limited flows [Costa and Macedonio 2005b;Pinkerton and Wilson 1994] will stop (i.e. flow length). Our results suggest that there is a highly non-linear temperature dependence, especially near the rheological cutoff where viscosity can increase by orders of magnitude over a narrow temperature range. Flow velocities can be estimated from a modified Jeffreys' equation [Jeffreys 1925;Takagi and Huppert 2007]: where v is flow velocity, ρ is density, g is gravity, α is ground slope, d is flow thickness, η is viscosity, and K is a geometric factor (K = 0.23 for a semi-circular channel geometry). Other channel geometries require numerical simulation beyond the scope of this study [e.g. Lev and James 2014]. Velocity profiles for Nyiragongo lavas were then calculated for ground slopes between 0-25°and apparent Newtonian viscosities from Table 2, for various flow thicknesses ( Figure 12). Dramatic rheological changes at our lower temperature and highest crystal fraction experiments suggest that extrapolations to temperatures, apparent viscosities, or crystal contents beyond our experimental conditions should be done with caution. Flow thickness has to be sufficient to overcome the apparent yield strength of the lava [Hulme 1974]. This relationship can be expressed by: where σ y is apparent yield strength and other symbols are as above. Figure 8 shows the different yield strength fits to our data compared to the predicted value from Dragoni [1989]. While the modeled values tend to underestimate yield strength, it does show that yield strengths are not expected to develop until much higher crystal fractions than we could achieve. The low crystal content of the experimental samples, particularly Nyiragongo lavas, result in negligible apparent yield strengths (Figure 8) during initial cooling, suggesting flows would be very thin. We note that Nyamuragira lavas have a higher σ y for a given temperature than Nyiragongo lavas. This is borne out in field observations of flow thicknesses. Nyiragongo lavas from the 2002 eruption range in thickness from 10 cm to~2 m [Allard et al. 2002;Baxter et al. 2002;Komorowski et al. 2002] while Nyamuragira lavas vary from 1-10 m aver-aging~3 m [Albino et al. 2015;Pouclet 1976]. These differences in σ y and flow thickness affect flow velocity in the opposite direction to the viscosity differences between the two lavas. The range of ground slopes roughly corresponds to the gradient along the main fissure system of the 2002 eruption of Nyiragongo. Figure 12 suggests that viscosities would have to be of the order of 10 1.3 -10 2.4 Pa s (~1220-1170°C) in order to reach velocities approaching 10 m s −1 for a 5°slope. Figure 12C and 12D provide a comparison of Nyamuragira lavas to those of Nyiragongo over the same range of slopes. For Nyamuragira lavas to reach a similar velocity (assuming a similar slope and thickness), viscosity would have to be of a similar order of magnitude (10 1.3 -10 2.4 Pa s), which would require much higher temperatures (~1260-1200°C). Bubbles may also influence the velocity by means of affecting the viscosity [Llewellin et al. 2002;Llewellin and Manga 2005] but this is not considered here. For another comparison, field estimates of channelized flow velocities from fissure 8 during the 2018 eruption of Kilauea ranged from 4 to 11 m s −1 [Dietterich et al. 2018], averaging about 8 m s −1 over the course of the eruption [Dietterich, pers. com.]. The lava flow velocity estimated for the 1977 eruption may drastically vary depending on the source of information. Pottier [1978] and Ueki [1983] give velocities of Figure 12: Comparison of flow velocities for Nyiragongo (top) and Nyamuragira (bottom) flows as a function of viscosity and slope. Velocities have been estimated using Jeffreys' equation (Equation 11). Newtonian viscosities used were taken from Table 2 and channel depths were estimated from field observations made during sample collection. Three different channel depths were used for each composition. NYI: 0.5 m (solid lines), 1 m (longdashed lines), and 2 m (short-dashed lines); NYA: 1 m (solid lines), 3 m (long-dashed lines), and 6 m (short-dashed lines). Colors represent the same viscosity for each channel depth with the solid lines labeled. The gray bar (top right) represents the estimated velocities from eye witness accounts of the Nyiragongo flows which are similar to the velocities reported for Hawaiian lavas.
roughly 4-6 m s −1 , based on local testimonies of the evolution of the eruption. Other publications provide slightly higher values of 5-11 m s −1 [Hamaguchi 1989;Hamaguchi et al. 1992] or 5-17 m s −1 [Durieux 2002b]. Tazieff [1977] provides a lava flow speed of 28 m s −1 , based on the fact that the lava flowed "at the speed of a car," but this estimate seems less realistic. For the 2002 eruption, the lava flows within channels may have reached speeds of a few m s −1 , with a mean speed of 0.5-1 m s −1 for the duration of the eruption . While probabilistic models are regularly used to produce useful hazard maps [e.g. Chirico et al. 2009;Favalli et al. 2009], the addition of rheological information could provide more robustness to this method.
Not only would this inclusion help to better constrain terminal flow lengths (i.e. cooling limited flows) but the timing of flows (i.e. flow velocity) can also be addressed. The temperature dependence of viscosity in particular can be very important as small changes in temperature can result in drastic changes in viscosity, especially at lower temperatures ( Figure 5). Other studies have attempted numerical simulations of lava flow paths at Nyiragongo Favalli et al. 2006;Favalli et al. 2009] only considering flows from the fracture system. However, a few cones exist in the northern part of Goma [Chakrabarti et al. 2009a] that are neglected as potential vent locations. These possible vents pose as great of a threat to the people of Goma Volcanica 3(1): 1 -28. doi: 10.30909/vol.03.01.0128 as lavas erupting from the fissures and should be considered in modeling.

Conclusions
Despite their broad chemical similarity, lavas from Nyiragongo and Nyamuragira have distinct rheological behaviors during cooling and crystallization. In this instance, the stratovolcano-Nyiragongouncharacteristically erupts lower viscosity lavas than the nearby shield volcano, Nyamuragira (presently). Nyiragongo lavas have a lower liquidus temperature (~1220°C) and begin crystallizing after a larger degree of undercooling (~75°C) than those of Nyamuragira. There is little crystallization that occurs over the range of experimental temperatures in either composition (NYI 1145°C, φ c = 0.04 ± 0.02; NYA 1151°C, φ c = 0.11 ± 0.02) suggesting that the increase in viscosity is mainly due to the temperature dependence rather than any physical or chemical effects of crystallization. No detectable yield strength was observed in either sample at crystal fractions below 10 vol. %, which agrees with many other studies suggesting yield strength development occurs at much higher crystal fractions than we were able to achieve. Nyiragongo produces particularly dangerous flows because of how fast the lavas can move even at large distances from the vent. Lavas remain fluid due to the minor crystallization occurring upon cooling below the liquidus. Viscosity and yield strength development remain relatively unchanged over the range of undercooling before >5 % crystallization occurs. The rheological transition between 1160°C and 1140°C (i.e. where the crystallization spikes in temperature space) means lava flow movement will be inhibited shortly after reaching this temperature range. Nyamuragira lava was 0.5-1.5 orders of magnitude more viscous than the Nyiragongo lavas over the experimental temperatures. Moreover, Nyamuragira lava starts to crystallize at higher temperatures at similar rates with respect to decreasing temperature (e.g. see similar slope on Figure 11B), whereby the rheological cutoff temperature is approximately 30°C higher for Nyamuragira lava compared to Nyiragongo lava (Figure 11). This allows Nyiragongo lavas to maintain higher flow velocities for longer periods of time than Nyamuragira lavas for a given eruption temperature and flow geometry. The comparison of the evolution of rheological properties of lavas erupted from both volcanoes demonstrate that Nyiragongo lavas possess high volcanic hazard potential. Thus far, volcanic hazard modeling has been insufficient in providing an all-encompassing picture of the threats posed by Nyiragongo to the city of Goma. There are few places in the world where a volcano poses such an imminent threat to a highly populated, urban area and this work attempts to better quantify the physical properties of lavas that may threaten the city.