The influence of emissivity on the thermo-rheological modeling of the channelized lava flows at Tolbachik volcano

.


Abstract.
The application of thermo-rheological models to forecast active lava flow emplacement and quantify important eruptive parameters of older flows has become a more common over the last decade.With the modification and adaption of these models to modular computing languages, they are now easier, quicker, and are being incorporated into studies of both terrestrial and planetary volcanism.These models rely on certain assumptions and input parameters, some of which such as emissivity are not well understood for molten materials.Without a well-grounded knowledge of how this parameter governs radiant cooling, remote measurements of temperature and models such as FLOWGO that rely on these temperatures to track cooling with time will be in error.Here, we perform a detailed FLOWGO-based modeling study of lava flows emplaced at Tolbachik volcano during the 2012-2013 and the 1975-1976 eruptions.Specifically, we have modified the FLOWGO model to incorporate a two-component emissivity model linked to the fraction of molten lava and cooled crust.We focus first on the large Leningradskoye Flow emplaced at the start of the 2012 eruption, relying on data from numerous other orbital sensors including MODIS, ASTER and ALI to constrain some of the model input parameters.The twocomponent emissivity adaption produced better fits to the final flow length, directly related to the crust cover percentage.We then applied the constrained model to the large Cone II flow formed in 1975, for which no satellite-based data are available.Results revealed that a nearly identical model fit was achieved with initial effusion rate of 700 m 3 /s or 1250 m 3 /s.However, for the higher the effusion rate, a lower the crust cover is needed to fit the flow length and width.This represents the first study to implement two-component emissivity into thermo-rheological modeling of lava flows.The results show that this is an important factor for model accuracy and

Introduction.
Basaltic volcanism is ubiquitous on Earth and the other inner solar system bodies.Over half of the world's volcanoes consist largely of basaltic-dominated systems occurring at every tectonic setting and on every continent [Walker, 2000].The basaltic volcanic eruptions at Tolbachik in Russia (2012Russia ( -2013)); Bardarbunga in Iceland (2014); Etna in Italy (2018); Piton de la Fournaise on Reunion Island (2018) and Kilauea in Hawaii ( 2018) reinforce the recurring hazard potential of basaltic activity.The 2018 activity at Kilauea formed a very large, stable channel from the vent to the ocean entry and destroyed numerous homes and property before coming to an end [HVO, 2018;Neal et al., 2018].Until the Kilauea eruption, the 2012-2013 eruption of Tolbachik volcano was the most thermally intense flow-forming eruption in the past 50 years, producing ~2.5 times more emitted energy than that of a typical eruption at Etna [e.g., MODVOLC, 2013;Pieri et al., 1990].Monitoring flow propagation direction, velocity and effusion rate, therefore, are critical for the flow models that have evolved over time.Several of these are focused on heat loss and down-flow topography to predict flow advance [e.g., Dragoni 1989;Favalli et al 2005;Garel et al 2014].In addition to topography, the dominant (internal) factors controlling flow propagation are the discharge rate combined with cooling and increasing viscosity [e.g., Walker 1973;Miyamoto and Papp 2004;Harris and Rowland 2009].All of these models, however, rely on surface temperature, a key source term parameter that is commonly measured using satellite-or ground-based thermal infrared (TIR) instruments [e.g., Flynn and Mouginis-Mark 1992;Wright and Flynn 2003;Donegan and Flynn 2004].It is the cooling of the flow's uppermost radiating glassy surface that is directly imaged by these TIR instruments.
Understanding the emissive properties of this surface thus becomes critical for any measurement or model reliant upon accurate knowledge of the kinetic temperature [e.g., Lillesand and Kiefer 1987;Ball and Pinkerton 2006;Harris 2013].
Thermo-rheological models of basaltic lava flows show that their morphological and dynamic evolution are governed by the interaction between the hot viscous core and the outer crust [e.g., Kilburn, 1993;Miyamoto and Sasaki, 1997;Miyamoto and Crown, 2006].These models were developed to examine flow evolution, heat loss, and ultimately their spreading rate, advance velocity, inundation area and flow front arrival time [Harris and Rowland, 2001;Keszthelyi et al., 2000;Vicari et al., 2007].A lava flow initially dissipates most of its heat radiatively, and with time and distance, forms a cooler glassy surface that will thicken, increase in viscosity, and eventually become a brittle crust [Hon et al. 1994].With a constant lava discharge rate, the increase in crust thickness and flow viscosity will eventually force a flow to stop (e.g., a cooling-limited flow) [Guest et al. 1987;Rhéty et al. 2017].Continued effusion upstream can then produce flow inflation and/or new break-outs and flow directions [Peterson et al. 1994;Crown and Baloga 1999].The magnitude of a flow's radiative cooling is related most strongly to the surface temperature and percentage of insulating crust [Flynn and Mouginis-Mark 1994], whereas the efficiency of that cooling is directly proportional to the emissivity of the hot fraction of the lava's surface [Holman 1992;Ramsey andHarris, 2013, 2016].
Emissivity is the unitless, wavelength-dependent fundamental property of a material and is sensitive to its composition, state and structure [e.g., Crisp et al. 1990;Kahle et al. 1995;Burgi et al. 2002] because it is directly related to the vibrational motion of the atomic bonds within the material (i.e., the petrology and structural state of the material).It is also influenced by the micron-scale surface roughness [Ramsey and Fink, 1999], and to a lesser degree, the look angle between the instrument and surface [Ball and Pinkerton 2006].Infrared spectra acquired remotely may be used to distinguish bulk wt.% SiO2, the presence of volcanic glass and the surface texture/vesicularity [e.g., Moxham 1971;Crisp et al. 1990;Ramsey et al., 2012].
Changes in material state (i.e., solid vs. molten vs. amorphous) or structure (i.e., composition) dramatically affect this infrared property.If emissivity is equal to unity at all wavelengths, the material is said to be a blackbody or perfect radiator.
The radiative temperatures of flows derived from thermal cameras or satellite instruments rely explicitly on knowledge of the surface emissivity, which is typically assumed close to unity [see review in Ramsey and Harris, 2013].Past anecdotal or poorly-constrained field measurements hinted at the fact that a melt's emissivity was lower than the cooled surface [e.g., Abtahi et al., 2002].If true, then so too is the derived effective radiation temperature [Pieri et al, 1990].An incorrect overestimate of emissivity will, therefore, overestimate the calculated radiative heat loss [e.g., Harris and Rowland, 2001;Keszthelyi and Denlinger, 1996;Keszthelyi and Self, 1998], which has a direct consequence on the cooling rate and the modeled final runout distance.Conversely, for multispectral TIR measurements, an incorrectly derived temperature results in an emissivity spectrum that is distorted in shape and incorrect in magnitude limiting accurate estimates of composition [Rose et al., 2014].
Here, we examine the effect of different emissivity values for the molten and crust fractions on the final length of channelized flows.More specifically, we use visible and infrared satellite data during the 2012-2013 eruption of Tolbachik, Kamchatka (Russia) to constrain certain modeling parameters used by PyFLOWGO [Chevrel et al., 2018].These satellite data provide important knowledge of the time-averaged discharge rate (TADR), channel width, radiant emission, and fraction of crust, all of which are used to refine the model results.We then use the constrained model to examine the large Cone II Flow that was emplaced in the same region 36 years earlier during the "The Great Tolbachik Fissure Eruption" (GTFE), which started on 6 Jul 1975 and ended on 10 Dec 1976 [Fedotov et al., 1991].

The Tolbachik eruptions.
The Tolbachik complex (55.83°N, 160.33°E) is located on the Kamchatka Peninsula, Russia and is comprised of two volcanoes: Plosky ("flat") Tolbachik (PT) with an elevation of 3,085 m and Ostry ("sharp") Tolbachik (OT) at 3,682 m.The westernmost OT is a peaked stratovolcano and the easternmost PT is a flat-topped shield volcano with a collapse caldera that formed during the GTFE eruption.To the south of these peaks lie ~875 km 2 of basalt flows, pyroclastic deposits, and a NNE alignment of cinder cones called Tolbachik Dol (TD) or "valley" [Fedotov et al., 1991].
The GTFE was Kamchatka's largest volumetric basaltic eruption in historic times.The axial portion of TD concentrates in a narrow (3 -4 km) zone where ~ 80% of all the eruptive centers are located (Figure 1).This eruption created a 20 km long chain of new cinder cones and flows, with the largest (the Cone II Flow) emanating from the second cone of the northern vents in the crater chain [Fedotov et al., 1991].This flow formed between 6 July and 10 Sept 1975, with its in the cracks on the crests of the large flow folds [Wessels et al., 2005].after the 50 th anniversary of the Institute of Volcanology and Seismology (IVS).Excellent details on the chronology, style and character of the eruption and eruptive products are given by Belousov et al. [2015] and Melnikov and Volynets [2015].The Kamchatka Volcanic Eruption Response Team (KVERT), which is part of IVS, initially reported episodes of volcanic tremor as early as 7 Nov 2012 [BGVN, 2012].That same day, observers from the village of Kozyrevsk (~ 40 km W) reported ash explosions and a red glow seen in the same area as the northern vents from the GTFE.Basaltic lava effused from two fissures along the west side of TD and a large thermal anomaly was immediately detected in polar orbiting, low spatial resolution satellite data [KVERT, 2012], which triggered later data acquisition from the high spatial resolution sensors used in this study.Vigorous fire fountaining activity produced channel-fed, fast-moving lava flows (Figure 2).By the time of the first cloud-free high spatial resolution orbital data acquisition by the Advanced Land Imager (ALI) on 1 Dec 2012, the primary flow (later named the Leningradskoye flow) emanating from the Naboko vent had already reached 11.3 km in length (Figure 3), eventually growing to ~ 14 km on 3 Dec, and ~ 17 km by 8 Dec with numerous smaller breakouts and new flows occurring over the subsequent months.The Leningradskoye flow is channelized, 15 m thick a'a flow [Belousov et al., 2015].It advanced at ~200m/h fed by the high effusion rates in the initial stage of the eruption.During this early phase of high effusion rates, the flow initially traveled through a deep and narrow channel with a velocity of 2-3 m/s [Belousov et al., 2015].
The flow within the channel was relatively uncrusted, which changed several weeks later as the discharge rate dropped and a short tube formed at the upper part of the channel near the vent.
Lava emerging from the tube at this point in time was nearly completely covered by a flexible, frothy crust 5-10 cm thick [Belousov et al., 2015].decreasing to 140 m 3 /s two weeks later [Gordeev et al., 2013;Dvigalo et al., 2014].This compares to the time-averaged discharge rate (TADR) derived from Moderate Resolution Imaging Spectroradiometer (MODIS) data of 278 m 3 /s on 1 Dec 2012, decreasing to an average of 102 m 3 /s in the following two weeks (Figure 4).It should be noted that the TADR values can have up to 30% error and represent a discharge rate averaged over the 12-24 hours prior to the image acquisition [Coppola et al., 2010].TADR were also calculated throughout the eruption following each cloud-free high resolution satellite image.Ramsey and Harris [2016] mapped the change in flow area using the high spatial resolution satellite data and a constant flow thickness of 5m (based on field reporting at the time) [e.g., Gordeev at al., 2013].The volumes are reported in Table 1.Throughout 2013, the flow field expanded later being redirected to the eastern side of the vent chain, eventually building three flow fields (Vodopadnoye to the NW, Leningradskoye to the SW, and Toludskoye to the SE) that covered over 35 km 2 [Dvigalo et al., 2013].1. Results of the first phase of the Tolbachik eruption using high spatial resolution satellite data to map flow field area and calculate effusion rate assuming a constant flow thickness of 5m.These results, calculated in real-time with each new satellite image acquired, compare quite well with later studies of the flow dimensions [Dvigalo et al., 2013;Gordeev et al., 2013].
Information on the erupted lava petrology is an important parameter constraint for the PyFLOWGO modeling.Petrologic information for the 1975 GTFE eruption comes from Fedotov et al. [1991] with the lava flow from Cone II falling into their "northern vents" category.These lavas are classified as magnesian basalt with moderate alkalinity, having a range of eruption temperature of 900 -1050 °C.Field-based estimates of lava viscosity made from lava flow front velocities ranged widely from 10 4 -10 10 Pa•s [Vende-Kirkov, 1978;Fedotov et al. 1991].One would expect a composition other than basalt for all but the very low end of this range, which was not the case.The northern vent lavas had an average bulk SiO2 wt.% of ~ 50 and 20 vol.
% phenocrysts, primarily clinopyroxene and olivine (3-8 mm).In contrast, Belousov et al. [2015] reported that the lavas from the 2012 TFE-50 eruption were more evolved with a slightly higher SiO2 content of 52 -53 wt.%.This, combined with the higher alkalinity, classifies these lavas as basaltic trachyandesites ( Reported LOI is assumed equal to H2O because of the basaltic rock composition.

High spatial resolution satellite data
The Advanced Land Imager (ALI) was launched aboard the Earth Observing-1 (EO-1) spacecraft in November 2000, and was deactivated by NASA in 2017.ALI had ten spectral channels in the visible/short wave infrared (VSWIR) region (0.43 -2.35 µm) with nine having a spatial resolution of 30 m and one (the panchromatic channel) spanning 0.48 -0.69 µm with 10 m resolution.The data were collected in four swaths, which spanned 37 km in width, with each swath having a length that varied from 42-185 km [Digenis et al., 1998;Hearn et al., 2001].ALI acquired numerous observations of the flows throughout the 2012-13 eruptive phase, including the first clear higher resolution image on 1 Dec 2012.
The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) was launched on the Terra spacecraft in December 1999 and measures ground-leaving radiance in 14 spectral channels [Yamaguchi et al., 1998].With the failure of the shortwave infrared (SWIR) subsystem in 2009, only two subsystems remain, each with different spatial/spectral resolutions: 299 the visible/near infrared (VNIR) sensor with three channels (0.56-0.81 µm) at a 15 m spatial 300 resolution and the TIR sensor with five channels (8.2-11.3µm) with 90 m spatial resolution.301 ASTER was tasked to acquire data of the Tolbachik eruption at every observational opportunity 302 (both day and night) from early Dec 2012 until late June 2013 (Figure 5).303 The ASTER science team also produced two versions of a global DEM (GDEM) dataset, which averages all the ASTER scenes acquired over a given region.This process greatly reduces instrument noise related errors and removes most cloud artifacts [Fujisada et al., 2005].16

High temporal resolution satellite data
The These are typically used to detect high-temperature volcanic thermal anomalies [Wright et al., 2002;Coppola et al., 2016].

MIROVA System
MIROVA (Middle Infrared Observation of Volcanic Activity) is an automated global hot spot detection system (http://www.mirovaweb.it)based on near-real time MODIS data [Coppola et al. 2016].The system completes automatic detection and location of high-temperature thermal anomalies through a series of processing steps.The system provides a quantification of the Volcanic Radiative Power (VRP in Watts), by means of the MIR method of Wooster et al. [2003]: where npix is the number of hot-spot contaminated pixels, Apixel is the pixel size (1 km 2 for MODIS) and RMIR,alert and RMIR,bk are the pixel-integrated MIR radiances (at 3.96 µm) of the i th alerted pixel(s) and background, respectively.According to Wooster et al. [2003], the constant of proportionality (18.9 m 2 ⋅µm⋅sr) in the equation allows estimations of VRP (± 30%) from hot surfaces having temperatures ranging from 600 -1500 K.This makes the application of eq. 1 to the lava flow surface a practical way to estimate the radiant flux of the active portion of the flow field, with essentially no contribution from the cooler and cooling flow areas.The cumulative Volcanic Radiative Energy (VRE) in Joules is calculated as the trapezoidal integration of the VRP time series.
During effusive eruptions, the calculation of Time Averaged Discharge Rate (TADR) of lava using satellite thermal data is an important input parameter for later modelling lava flow advances and dynamics [Harris et al., 2016].This approach, commonly referred to as the "thermal proxy", is derived from the contribution of several studies that have continuously revised and refined the theoretical framework and its practical application to real cases [e.g., Pieri and Baloga 1986;Crisp and Baloga 1990;Harris et al., 1998;Wright et al., 2001;Harris et al., 2007a;Harris and Baloga 2009;Dragoni andTallarico 2009, Coppola et al, 2013;Garel et al., 2012, 2014, 2015, Tarquini et al., 2017].A specific analysis of the TADR calculations using MODIS-MIROVA data [Coppola et al., this volume] outlines the methods and limits of this approach during the emplacement of large lava flows, such as the case for the TFE-50 Tolbachik eruption.Following Coppola et al. [2013], a single coefficient, called radiant density (crad in J m -3 ), can be calculated to describe the appropriate relationship between radiated energy (VRE) and erupted volume (Vol) for the observed eruption: For the Tolbachik eruption we considered a total volume of the lava flow of ~573 x 10 6 m 3 [Dai and Howatt, 2017] and calculated a crad equal to 1.08 × 10 8 J m -3 , a value typical of lava flows having mafic composition [Coppola et al., 2013].It follows that for any VRP measurements, the corresponding TADR can be calculated following: Sensor webs The improved temporal coverage of the Tolbachik TFE-50 eruption by the ALI and ASTER sensors was made possible by detection and communication protocols, as well as software to implement those protocols, all designed to create "sensor web" networks [Davies et al., 2016;Ramsey et al., 2016].Both ASTER and ALI have 16-day nominal repeat observational frequencies at the equator, which improves at higher latitudes due to the convergence of overlapping orbit swaths.However, with the sensor web programs, these time periods are greatly improved [see Davies et al., 2016;Ramsey et al., 2016].
The sensors on the EO-1 spacecraft were used as part of the Volcano Sensor Web (VSW), which relied on onboard software for rapid scheduling and data processing using reports of volcanic unrest as triggers [Davies et al., 2016].The VSW did not rely solely on detections from other orbital sensors, but also ground-based detection and communicated reports.The rapid retasking of the spacecraft, off-nadir pointing as well as the use of onboard data processing and downloading commonly enabled very rapid response times to volcanic eruptions.
A program that links high temporal/low spatial resolution data detection of thermally-elevated anomalies to specific scheduling of ASTER data has been in place since 2004 [Ramsey and Dehn, 2004;Carter et al., 2007Carter et al., , 2008;;Rose and Ramsey, 2009;Duda et al., 2009;Ramsey, 2016].The urgent request protocol (URP) global system started with (and has been heavily focused on) the northern Pacific region [Duda et al., 2009;Ramsey and Harris, 2013;Ramsey, 2015].The URP program allows data acquisition frequency as high as night-day-night observational triplets acquired approximately every two weeks and day-night pairs acquired approximately every five days during the Tolbachik eruption.
The presence of the VSW and the URP provided a much larger volume of higher spatial resolution data than would normally be available for such an eruption.For example, ALI acquired the first high resolution image of the Tolbachik eruption on 1 Dec 2012, whereas the first clear ASTER observation came ~ 30 hours later during a nighttime overpass.ASTER was also able to acquire a daytime scene only 13.5 hours later due to the orbit configuration of the Terra satellite [Ramsey, 2016].During the first seven months of the Tolbachik eruption, there were 48 cloud free/nearly cloud free acquisitions by all high resolution satellite instruments (an average of one scene approximately every 4 days).Of the total number of scenes, 33 were ASTER, 11 were ALI and 4 were from other sensors.The data acquired in the first few weeks of the eruption in particular provided some of the best synoptic views, documenting vent locations, flow lengths/direction, channel widths and their changes over time (Figure 5).

PyFLOWGO modeling
The original one-dimensional FLOWGO model tracks the thermal and rheological evolution of a lava control volume using a series of heat loss equations coupled with the Jeffreys equation for Newtonian flow in an open channel, modified for a Bingham fluid [Harris and Rowland, 2001;Harris and Rowland, 2014].Lava is tracked down-flow with its cooling, crystallization, viscosity and yield strength re-calculated at each down-flow step to estimate the velocity [Harris et al., 2007b;Riker et al., 2009;Robert et al., 2014].The model was applied to various cooling limited lava flows [Harris and Rowland, 2001, Harris et al., 2005, Harris and Rowland, 2012, Rhéty et al., 2017].PyFLOWGO [Chevrel et al. 2018] implements the model using the Python open source software [ van Rossum, 1995].
To correctly simulate the evolution of down-flow lava properties, the user is allowed to change the input parameters as well as the thermo-rheological models (i.e.heat loss mechanisms, crystallization rate, temperature-and crystallinity-dependent viscosity, crust cover fraction, etc.) within plausible limits [e.g.Harris et al., 2007b] so as to best-fit a natural flow.Input parameters cover broad categories of channel dimension, eruption conditions, the radiative, conductive, convective parameters, and material properties of the lava.Within these categories are some terms (typically assumed) that are specifically related to (and derived from) the TIR radiance, including emissivity, lava/crust temperatures, and vesicularity.Typically, lava emissivity and vesicularity are assumed and held constant throughout a model run as the temperature of the crust cools.As originally specified by Harris and Rowland [2001], the fraction of this crust increases down-flow as flow velocity decreases (rather than being a function of the radiative cooling), implying a flow regime in which the crust is more stable at lower flow velocities (and presumably lower temperatures).
For thermo-rheological models such as PyFLOWGO that rely on radiant cooling, the question arises: to what degree does an incorrect emissivity (and perhaps other) assumption(s) affect the final model results.To first test whether emissivity plays a significant role in the radiant cooling, and therefore in determining the final cooling-limited flow length, we modified PyFLOWGO with a radiative heat flux (  ) model that uses an "effective" emissivity: Here, it is not just the effective radiation temperature (Teff) that depends on the crustal fraction (fcrust), but also the effective emissivity (  ).This allows the effective emissivity to be computed via a two-component emissivity model where the cooler lava crust emissivity (  ) is different than that of the molten lava emissivity ( ℎ ): This differs from previous FLOWGO modeling that assumed only a single emissivity ( ℎ =   = 0.95).Based on recent work by Ramsey and Harris [2016], Lee and Ramsey [2016] and Lee et al. [2010], the emissivity of molten lava in situ ( ℎ ) and silicate glasses in the laboratory is likely much lower (avg.~ 0.6) over the 5 -25 μm TIR region.For the emissivity of the crusted surface (  ), we use the same value (0.95) from prior studies [Harris, 2013].
Assuming this two-component emissivity model, therefore, results in a lower radiant heat loss for situations where the fcrust < 1.This should result produce a slower cooling flow and therefore, a greater cooling limited distance (maximum length).The lower emissivity of the molten fraction reduces its radiative efficiency; and hence, the control volume cools more slowly.

The two-emissivity model
To estimate the maximum difference in the cooling limited distance between the two emissivity models, we first ran PyFLOWGO with a constant slope of 5°, crust temperature and crust fraction down the entire flow length.The simulations are shown between an uncrusted flow (fcrust = 0) with an effective emissivity of 0.6 and a fully crusted flow (fcrust = 1) with an effective emissivity of 0.95; as well as a flow that is half crusted (fcrust = 0.5) having an effective emissivity of 0.78.The maximum difference between the two models is expected to be between the uncrusted and the fully crusted flows.Results show a linear relationship between crust cover/emissivity to the final modeled flow distance.The two-component emissivity model produces flows that are a maximum of 34% longer than flows using a single constant emissivity (Figure 7).The difference in cooling limited distance here is only dependent on the crustal fraction, which is a function of flow velocity and directly related (i.e., the complement) to the fraction of molten material having the lower emissivity.

Petrology
We analyzed the whole rock composition of two samples which were sampled from the lower, folded region of the Cone II flow and near the terminus of the Leningradskoye Flow (Table 2).
The Leningradskoye flow sample is broadly consistent with reported values [e.g., Gordeev et al., 2013;Belousov et al., 2015;Plechov et al., 2015].The results from the Cone II flow sample, however, do vary somewhat from those of the later stage northern vents reported by Fedotov et al., [1991].Our results show lower Al2O3 and higher in CaO and MgO, each by ~2 wt.%.These variations would still classify the sample as basaltic, but are more consistent with the Fedotov et al. [1991] results for the earlier stages of this eruption, which were emplaced before the Cone II flow.These compositions were used for the petrologically-dependant fluid viscosity model used by PyFLOWGO (Tables 2 and 3).

The 2012 Leningradskoye Flow (TFE-50 eruption)
Using MIROVA, we calculate that the TFE-50 eruption radiated approximately 5.87 (±1.76) × 10 16 J into the atmosphere during the course of the eruption.Application of eq. 3 to the Tolbachik MODIS data, suggests that the TADR peaked at the very beginning of the eruption reaching 300 (±100) m 3 s -1 on 29 November 2012 at 01:09 UTC.Later, the eruption showed a progressive decline in intensity, with effusion rates dropping to 278 m 3 /s on 1 Dec 2012, less than ∼100 m 3 s -1 on 13 December 2012, ∼50 m 3 s -1 , on 22 January 2013, and ∼5 m 3 s -1 by 22 August 2013 (Figure 4).These values correspond closely to those derived from visual mapping of the flow field using the higher resolution ASTER/ALI data (Table 1).

MODEL NAME MODEL CHOICE REFERENCE
"crystallization_rate_model" basic Harris and Rowland (2001)   m/pixel (ASTER VNIR) and, second, to the intense radiance, which where convolved with an instrument's point spread function, causes excess radiance to be detected in the surrounding lava-free pixels.This results in the channel-related thermal anomaly to appear wider than the actual widths measured in the field.The effect was most noticeable in the upper part of the flow where the emitted radiance was most intense.The channel width at the vent was therefore set to 30 m, the minimum resolvable distance in the ALI data and the ASTER DEM.To match the effective MIROVA effusion rate over the at-vent slope, the depth is set at 6.1 m, which is in agreement with field-based estimates [Belousov et al., 2015].Some PyFLOWGO input parameters were obtained by information obtained during the eruption [e.g., Plechov et al., 2015, Volynet et al., 2015;Gordeev et al., 2015;Belousov and Belousova, 2018] or by assumptions based on other well-constrained basaltic flows following Harris and Rowland [2001] and Chevrel et al. [2018].For example, to estimate the initial viscosity we consider the petrologic observations from Plechov et al. [2015], which include the interstitial glass composition having 0.03 wt.% H2O, and the cooled lava having an average of 25 vol.% crystals and 6 vol.% bubbles.We used the model of Giordano et al. [2008] for the interstitial melt viscosity in association with the Einstein-Roscoe model for computing the effect of crystals.Because the vesicle fraction is as low as 6 vol.%, the effect of bubble on viscosity is neglected here.Model results produced an initial viscosity of 1.9 × 10 4 Pa•s at the eruption temperature of 1082 °C, which is in agreement with field estimations and measurements done by Gordeev et al., 2015 andBelousov andBelousova, 2018].The simulation that best fit the flow length constrained by the channel widths and emitted radiance down-flow measured from the ALI image was produced using the parameters reported in Table 3 and is presented in Figure 8.The best fit model result used a variable crustal coverage as part of the "basic" model in 629 PyFLOWGO, originally proposed by Harris and Rowland [2001].This model allows   to 630 increase as a function of mean flow velocity (  ): 631 where   , the initial (at vent) crust fraction and α, a coefficient, were derived for a poorly 633 insulated flow by Harris and Rowland [2001] to be 0.9 and -0.16, respectively.The velocity of 634 the lava within the channel varied between 2 and 3.5 m/s within the first kilometer and then 635 progressively decreased down-flow (Figure 9d).These values are also in agreement with fieldbased estimates reported in Belousov et al. [2015].In these simulations, the difference between the single and the two-component emissivity is small because the modeled crust coverage fraction is high, varying from 60% -90%.(Figure 10).This range is somewhat higher than seen in the field images of the open channel, which varied between 35% -60% (Figure 2).However, the field-based values are derived from the visibly darker fraction of crust on the channel surface.It is possible that the "radiative crust" (i.e., surfaces that are still visibly red but have cooled enough to raise the emissivity) is higher.(Figure 6a) from the ALI band 7 (0.87 668 μm).The subpixel hot fraction emitted 669 radiance was estimated using the two-670 component mixing approach of Dozier 671 [1981] and then compared to the 672 PyFLOWGO model output (Figure 8b).673 The results show that the two-component emissivity model does result in a spectral radiance 674 that is closer to that of the ALI data.However, the fit deteriorates down-flow with the measured 675 radiance decreasing to a near constant value, approximately an order of magnitude lower than 676 the PyFLOWGO model predicted values.This overestimation by the model could be due to 677 several factors.The ALI data are resampled during the conversion from the original radiance 678 data (L1R) to the distributed L1Gst data product.This could cause mixing of radiant energy from 679 adjacent pixels containing cooled lava or snow, thus lowering the pixel-integrated values.Furthermore, the ALI data have not been atmospherically corrected.There were significant steam plumes and clouds encroaching on the upper part of the flow and smaller vapor plumes are seen emanating from the lava channel along the entire length.This atmospheric interference will lower the pixel integrated radiance [see Sawyer and Burton, 2006].Finally, whereas PyFLOWGO calculates emitted radiance based on the hottest lava temperature at a particular point along the flow, the satellite data captures emitted energy from all temperatures within a given pixel.It is very likely that within a 30 m pixel, there exists some fraction of hot material plus that of cool/cold material, as well as numerous surfaces at temperatures between these two endmember temperatures [Dozier, 1981;Rothery et al., 1988;Harris, 2013].These would not be resolved in the two-component subpixel analysis and yet would lower the effective pixelintegrated spectral radiance.Therefore, use of spectral radiance as a constraint on PyFLOWGO model output is best applied to higher-resolution (i.e., ground-or airborne) data together with in situ measurements of the intervening atmosphere.This flow extended 5 km and was characterized by a very wide channel, oversteepened levees and a folded terminus with a 50 m high flow front.The morphology of the flow makes it appear perhaps more silicic than the actual composition, which is magnesian, moderately alkaline basalt [Fedotov et al., 1991].Samples acquired from this flow in 2005 were found to be basalt with 50.5 wt.% SiO2 (Table 2).Although the composition is basaltic, the effective viscosity was estimated in the field to range between 10 4 -10 7 Pa•s and up to 10 9 Pa•s, based on lava flow front velocities [Vende-Kirkov, 1978;Fedotov et al. 1991].
To simulate the evolution of the thermo-rheological properties down-flow, we used the same input parameters as used for the 2012 Leningradskoye lava flow but imposed an initial melt viscosity on the low-end (5 x 10 4 Pa•s) of the range provided by Vende-Kirkov [1978] together with an eruption temperature of 1050 °C [Fedotov et al., 1991].We were able to fit the channel width and the final length with an effusion rate of 700 m 3 /s and an initial channel dimension of 70 × 8.7 m (Figure 11).The difference between the two-component emissivity model and the single emissivity model was negligible considering the high initial crust cover fraction of   = 0.9 used (Figure 12).However, the large channel and thickness of this flow indicates a higher initial effusion rate and flow velocity could have been possible.In a case like this, a similar model  The PyFLOWGO model provides a great degree of flexibility on many of the input variables and choice of viscosity models.This flexibility also constraining the model output quite difficult without knowledge of at least some of the input parameters, and highlights the need for independent flow parameter data against which to best-fit model-driven scenarios.Past studies using FLOWGO modeling commonly compared the results to well-constrained field data, such as the channel dimensions, measurements of the crust percentage, velocity and/or cooling profiles down-flow [Harris et al., 2007b;Wright et al., 2008, Harris and Rowland, 2009, 2012].
Here, we apply PyFLOWGO, constrained primarily by high spatial (e.g., ASTER and ALI) and high temporal (e.g., MODIS) resolution satellite data.The MODIS time-averaged discharge rate and the ASTER-derived DEM data proved to be the most important as they provide the volume flux source term, and topographic underlay required by the model, respectively.The ASTER/ALI data also allowed the flow length, channel widths, advance rate, temperature and spectral radiance to be constrained to some degree.Two parameters measured from the high spatial resolution data (channel width and spectral radiance) proved more difficult to constrain.The excessive radiance from the wider, less crustcovered channels closer to the vent resulted in overestimations of the measured channel widths.
Conversely, channelized regions with a modeled crust percentage greater than ~75% had emitted spectral radiance much lower than the PyFLOWGO predictions.Incomplete atmospheric correction, image resampling and pixel mixing likely produces emitted radiance values at scales too complex to be extracted from the 30 m pixels examined here; an issue that warrants future study.
Despite the paucity of constrained ground-based data and the complexities regarding higherlevel modeling of some of the high resolution image data, we were able to model the Leningradskoye lava flow by fitting the final flow length and starting discharge rate.Model output of core temperature, crystal fraction, mean velocity and viscosity are in good agreement with ground and satellite control.Applying these constraints then to the Cone II flow combined with data extracted from the ASTER GDEM allowed us to estimate the discharge rate as well as the effect of variable emissivity and initial crust cover.The difference between two-component and the single emissivity models was negligible where considering a high initial crust cover fraction, but significant for higher effusion rate with a lower initial crust cover fraction .
Such a validation approach is important for interpreting older flows where information on their emplacement is either limited or nonexistent.
The primary focus of this study was to apply PyFLOWGO using data constrained as much as possible by only satellite data, and in the process, test one variable (emissivity), which has always been assumed constant in past studies.We modified PyFLOWGO to accept two values for emissivity and assigned those to the crust and molten fractions in the channel.The effective emissivity term was thus a weighted sum of the crust fraction having a high emissivity (0.95) plus the molten fraction having the lower emissivity (0.60).This change has a measureable (linear) impact on the PyFLOWGO results for a theoretical flow and where applied to the actual Tolbachik flows resulted in a slightly improved fit for the 2012 Leningradskoye flow and significantly improved fit for the 1975 Cone II flow.As one might expect, the largest change in the model results arises where a flow initially enters the channel free or nearly free of cooling, higher emissivity crust.If a flow is modeled with a high initial crust fraction, the two end-member emissivity modification will only provide a marginally better final fit.However, what we did not explore here is assumption in PyFLOWGO that crust formation is related solely to flow velocity rather than temperature.A more complex solution incorporating crust formation with flow cooling (modulated by the lower emissivity) is warranted.Field images of the 2012 flow channels commonly showed significantly less crust at a given position along the flow than was predicted by PyFLOWGO.If this is the case, then we predict the lower emissivity of the molten fraction of the channel will have more of an impact.
The emissivity measurement of molten material has been presented in numerous recent studies both in the laboratory and the field [Ramsey and Harris, 2016;Lee and Ramsey, 2016;Lee et al., 2010;Abtahi et al., 2002].The assessment of Burgi et al. [2002] for the emissivity of the surface of the active lake at Erta Ale (0.74), albeit in the 1.1-1.7 µm region, is in agreement with this study [see also Flynn et al. 2001].This fundamental change in our understanding of how efficiently materials radiate heat prior to solidifying is important for models such as PyFLOWGO but also in any situation where an accurate thermal infrared non-contact temperature is made.Emissivity measured over multiple wavelengths produces a spectrum that relates information about the material's petrology, its state of atomic bond vibrations and angles, and the amount of glass formation [see Ramsey and Fink, 1999].The broadband emissivity, an average over the entire wavelength region, is a measure of how efficiently that material emits heat.It should be noted that the prior studied of molten emissivity were conducted in the TIR region (~ 5 -25 μm in the laboratory and ~ 8 -12 μm in the field/from orbit).Radiant heat loss, however, is occurring over a larger portion of the electromagnetic spectrum, most notably in the short-wave infrared region (~ 1 -5 μm).If the emissivity is closer to unity at these wavelengths, then the impact of the lower emissivity in the thermal infrared region will be mitigated to a degree.
Further laboratory studies of emissivity in the 1-5 and 8-14 μm regions are planned to clarify this question.Further study into the model's logic of linking crust formation to velocity rather than temperature is also planned.This is likely resulting in an underestimate of the effect of the two end-member emissivity approach.More model simulations with this logic changed to temperature or temperature/velocity dependency are planned.
central channel later covered by a smaller and thinner flow emplaced between 11 and 15 Sept 1975.The ~ 5 km long Cone II flow is 50 -60 m thick with well-developed 30-m-high levees and a central channel that disappears down-flow where it spreads over flatter topography.In this zone, dominant flow folds (15 -20 m high) are present perpendicular to the flow direction.Also present are regions of elevated thermal output that persist today.The slow cooling of this thick, voluminous unit produces surface temperatures in excess of 100° C, which have been measured

Figure 1 .
Figure 1.ASTER VNIR image acquired on 2 March 2013 centered on the Tolbachik Dol ("valley"), with channels 3 (0.807 μm), 2 (0.661 μm), 1 (0.556 μm) in red, green, blue, respectively.Spatial resolution is 15 m/pixel.The longest flow emplaced early in the eruption sequence is the Leningradskoye flow.Similarly, the other large flow analyzed here is the 1975 Cone II flow.Note the snow-free regions toward the flow's terminus due to elevated heat flow over 40 years after emplacement.Inset image is the TFE-50 vent region (brightened with a histogram equalization stretch) shown at full resolution.The active Naboko vent and channels shown in red are denoted by arrows.

Figure
Figure 2. Aerial photographs of photograph taken on the same

Figure
Figure 5.Time series of ASTER TIR data acquired over the first 3 weeks of the TFE-50 eruption.All images are shown at the same scale, with a guassian stretch applied.Spatial resolution is 90 m/pixel.The (n) or (d) after the date refers to data acquired at night or day local time, respectively.(A) 2 Dec 2012 (n).(B) 3 Dec 2012 (d).(C) 11 Dec 2012 (n).(D) 12 Dec 2012 (d).(E) 18 Dec 2012 (n).(F) 20 Dec 2012 (n).Images A -C and F contain various degrees of intervening cloud cover, though not enough to entirely block the TIR radiance, but significant enough to blur details and make accurate temperature retrievals impossible.Image D contains an inset of the VNIR data from the same date, which is able to resolve the individual open channels (shown by yellow arrows).
Version 1 (v1), released in 2009, was compiled from over 1.2 million individual scenes.The improved GDEM version 2 (v2) was released in 2011 adding 260,000 additional DEMs acquired from 2008 -2011, improving coverage and further reducing data artifacts.The refined production algorithm for v2 also provides finer horizontal resolution and increased accuracy (e.g., vertical accuracy of ~17 m on average)[Tachikawa et al., 2011].The ASTER GDEM v2 is used as the topographic base layer of choice for the analysis of the Cone II flow from the 1975-76 eruption.The creation of GDEM v2 ended with data acquired prior to the start of the TFE-50 eruption.Therefore, we rely on the single-scene DEM's for the 2012-13 eruption.Pre-flow topography is extracted by mapping the position of the 2012 flow channels from the ALI SWIR data and ASRER TIR data on the GDEM.These elevation data are used in PyFLOWGO (Figure6).Individual ASTER scene DEMs collected during the eruption were analyzed to determine channel width, flow thickness, and channel topography where possible.However, the presence of volcanic fume and cloud combined with the scene-dependent noise limited the use of the individual scene DEMs.Pre-flow topography for the 1975 flow was extracting from the GDEM using transects next to the flow and directly down the centerline of the main channel.Channel widths of this larger flow are also easily resolved in the GDEM data.Five channel cross sections were extracted down the length of the flow and used to constrain PyFLOWGO model runs.

Figure 6 .
Figure 6.High resolution orbital data draped over the ASTER-derived DEM produced from the VNIR image acquired on 28 Feb 2013.DEM posting = 30 m. (A) ALI image acquired on 1 Dec2012 with channels 9, 8, 7 in red, green, blue, respectively (see Figure3).White triangles denote extracted radiance locations (see Figure8B).(B) ASTER image acquired on 28 Feb 2013 with channels 3, 2, 1 in red, green, and blue, respectively.The active channel in the ALI image was traced, exported as a vector and shown in green.

Figure 7 .
Figure 7. Percent difference between the run-out distances of the PyFLOWGO simulations using a 2component emissivity model (ε hot = 0.6 and ε crust = 0.95) versus a single emissivity model (ε hot = ε crust = 0.95) as function of crust cover fraction.There is no difference in run-out distances if the flow is fully crusted.By contrast, if initially uncrusted, the maximum difference is 33%.Crust temperature plays a minor role, but those changes do not exceed the symbol size on the figure.

25
Because the ASTER GDEM is created from all ASTER scenes spanning from 2000 -2011, it does not resolve the TFE-50 eruption lava flow fields.The ASTER single scene DEM created from the 11 Jan 2018 acquisition, however, was relatively free of these data errors and represents the surface five years following the TFE-50 eruption.This DEM was precisely geolocated to the ASTER GDEM, which was then subtracted from it.The outlines of the flow fields are resolved as a positive elevation anomaly.Transects taken perpendicular to the Leningradskoye Flow show steep levees and flow fronts with maximum thicknesses of 35 m.In some of these transects, the channel is visible and ranges from 30-50 m in width and 5 -20 m in depth.It was not possible, however, to discern the channel at the vent, perhaps because it had been filled in by subsequent flow activity.A transect down this flow varies in average thickness from 46 m near the vent to 16 m over the rest of the flow.PyFLOWGO was run to simulate the Leningradskoye lava flow that formed between 29 Nov 2012 and 1 Dec 2012.In that period, the flow was fed by an open channel and reached a length of 11.3 km, which is 70% of its final length.After 1 Dec 2012, the effusion rate dropped and later lava tubes formed allowing the flow to extend to its final length (16.4 km).We did not attempt to fit the final flow length because PyFLOWGO is designed only for open channel systems.We use the effusion rate for this time period determined by MIROVA (278 m 3 /s) for model initialization and the steepest line of descent measured on the pre-flow ASTER GDEM.The active channel widths are measured directly from the ALI and ASTER infrared data.These are likely overestimates, however, due first to the minimum spatial resolution of 30 m/pixel (ALI) and 15 Figure 8. Best fit FLOWGO simulation for the Leningradskoye lava flow of the 2012 Tolbachik eruption using input models and parameters given in Table 3.The continuous line is for a two-component emissivity model, whereas the dashed line represents a singular emissivity of 0.95.The red dots are channel width measurements from ASTER and ALI images.The vertical red lines represent the flow length on 1 Dec 2012 measured from the ALI image.(A) Modeled channel width.Misfits near the beginning of the flow are likely due to the intense radiance from the lava in these regions causing the channels to appear larger than their actual width due to radiance bleeding into neighboring pixels.(B) Modeled emitted radiance.Misfits at greater distances here are likely due to increasing percentages of cooled crust/older lava mixing within these pixels.

Figure 9 .
Figure 9. Thermo-rheological variations for the 2012 Leningradskoye lava flow properties using the best fit simulation and the input models and parameters given in Table 3.The continuous line is the two-component emissivity model results, whereas the dashed line represents a singular emissivity of 0.95.The vertical red lines represent the flow length on 1 Dec 2012 measured from the ALI image.The singular emissivity produces model results that consistently over-predict the actual flow length.

Figure 10 .
Figure 10.Variations in the crustal coverage, The 1975Cone II Flow (GTFE eruption)     The Cone II flow dimensions are easily resolved in the ASTER GDEM making it much more straight forward to measure the channel width and depth down-flow.Six perpendicular transects were taken from the base of Cone II the point where the flow spreads laterally and folding masks the end of the channel.The channel averaged 140 m in width and 11.8 m in depth over the six transects.At the base of Cone II, the channel is between 60-90 m wide and 8-9 m deep.These values are used as bounding constraints for initiation of PyFLOWGO (

Figure
Figure 12.Variations in the crustal

Figure 13 .
Figure 13.Thermo-rheological variations for the 1975 Cone II lava flow properties using the best fit simulation.The input models and parameters are similar to those used for the 2012 eruption, except the few presented in Table 4.The green lines show the results for an effusion rate of 1250 m 3 /s and an initial crust fraction (  ) of 0.4, whereas the black lines represent an effusion rate of 700 m 3 /s and   = 0.9.The solid lines represent the two-component emissivity model results, whereas the dashed lines are a singular emissivity of 0.95.The vertical red lines represent the final flow length on 1 Dec 2012 measured from the ALI image.

Table 2
).Over time, the erupted lavas became more basic with SiO2 weight percentages decreasing.Plecheov et al. [2015]reports the average crystal content for the Naboko vent lavas, which produced the Leningradskoye flow, ranged from 23 -29 vol.% with a porosity of ~ 6 vol.%.The primary crystals were plagioclase, olivine and titanomagnetite.Based on geothermometry analyses, they estimate the eruptive lava temperature was ~1080 °C.

Table 2
Belousov et al. [2015]nalysis of major oxides done for this study, comparing one sample from the 1975 GTFE to one from the 2012 TFE-50 Tolbachik eruptions.The 2012 sample is higher in AlO3 as reported byGordeev et al. [2013]andBelousov et al. [2015], plotting in the same geochemical space as their sample analyses.Comparing the SiO2 wt.% to the sum of the K2O + NaO2 wt.%, classifies the 1975 sample as a basalt and the 2012 sample as a basaltic trachyandesite, consistent with prior studies.
MODerate resolution Imaging Spectroradiometer (MODIS) is carried on both the Terra and Aqua NASA-EOS satellites, which have been flown in sun-synchronous polar orbits since February 2000 and May 2002, respectively.The two MODIS instruments provide radiometric data in 36 spectral bands from 0.4 µm to 14.4 µm.Of these, 29 channels collect data in the IR region with a nominal spatial resolution of 1 km at nadir.The ±55-degree scan angle produces 2330 km swath widths with data acquired approximately four times per day for a given target.At the higher latitudes of the Kamchatka peninsula, however, 6-10 overpasses per day are common due to pole-ward orbit convergence.MODIS IR channels, include a middle infrared (MIR) and a TIR channel, centered at 3.96 µm (channels 21 & 22), and 12.02 µm (channel 32), respectively.

Table 3 .
All PyFLOWGO input models and parameters that were used for the Leningradskoye lava flow of the 2012

Table 4 .
Table 4).Specific PyFLOWGO input models and parameters that differ from those in Table 3 that were used for the Cone II lava flow of the 1975 Tolbachik eruption.