The influence of variable emissivity on lava flow propagation modeling

Modeling lava flow propagation is important to determine potential hazards to local populations. Thermo-rheological models such as PyFLOWGO track downflow cooling and rheological responses for open-channel, cooling-limited flows. The dominant radiative cooling component is governed partly by the lava emissivity, which is a material property that governs the radiative efficiency. Emissivity is commonly treated as a constant in cooling models, but is shown here to vary with temperature. To establish the effect of temperature on emissivity, high spatiotemporal, multispectral thermal infrared data were acquired of a small flow emplaced from a tumulus. An inverse correlation between temperature and emissivity was found, which was then integrated into the PyFLOWGO model. Incorporating a temperature-dependent emissivity term results in a ∼5% increase in flow length and < 75% lower total cumulative heat flux for the small flow. To evaluate the scalability of this relationship, we applied the modified PyFLOWGO model to simulations of the 2018 Lower East Rift Zone fissure 8 flow, emplaced between May 27 and June 3. Our model improves the emplacement match because of the ~ 30% lower heat flux resulting in a ∼7% longer flow compared to modeling using a constant emissivity (0.95). This 5–7% increase in length prior to ocean entry, realized by an accurate temperature-dependent emissivity term, is critical for developing the most accurate model of future flow hazard assessments, particularly if population centers lie in the flow’s path.


Introduction
Lava flow modeling is a powerful tool for quantitatively forecasting lava propagation and subsequently improving the accuracy and reliability of lava hazard assessments (Ramsey and Harris 2013). This importance is reinforced by recent eruptions at Kīlauea (Hawai'i), Piton de la Fournaise (La Réunion), Etna (Italy), and Pacaya (Guatemala), which produced lava flows that posed serious risks to local societies. For example, in the summer of 2018, activity at Kīlauea volcano emplaced numerous channelized lava flows up to 15 km long, destroying over 700 buildings in the Puna district (Neal et al. 2019). Flow propagation modeled during an ongoing eruption provides useful guidance to hazard response coordinators in populated areas so more informed risk reduction measures may be enacted .
Lava flow propagation models require at least some input parameters. The length and dispersion of a lava flow are controlled by the temperature, effusion rate, crystal/vesicle content, and the topographic slope. Temperature is the foremost variable because it influences rheology (Harris et al. 1998;Cashman et al. 1999;Gregg and Fink 2000). Indeed, the majority of basaltic lava flow propagation models assume that temperature is inversely related to the viscosity (Park and Iversen 1984;Dragoni and Tallarico 1994;Harris and Rowland 2001). In a hazard response scenario, gathering the requisite petrologic (if no prior data are available) to perform accurate flow modeling is time-consuming, and can take several days to process (Ramsey and Harris 2013). Sometimes, thermal properties are estimated or assumed from previous field or laboratory studies (e.g., Harris and Rowland 2001;Avolio et al. 2006;Bilotta et al. 2012), whereas other recent studies integrate satellite-based Advanced Very High-Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS) thermal infrared (TIR) data into flow modeling (e.g., Negro et al. 2008;Herault et al. 2009;Vicari et al. 2009Vicari et al. , 2011. Generally assumed, however, is the emissivity of the lava needed to derive the radiant temperature, heat flux, and cooling of the flow over time. Emissivity is a wavelength-dependent, unitless property that determines the efficiency at which radiant energy is emitted by a material. Most natural surfaces are not perfect emitters (i.e., they are not blackbody surfaces). Emissivity is controlled by the atomic vibrational motion of the material and the physical properties of the surface such as roughness and particle size. Therefore, it varies with composition, vesicularity, and surface expression (e.g., Christensen et al. 2000;Abtahi et al., 2002;Vaughan et al. 2005). However, temperature was never considered a controlling variable. Prior emissivity measurements quantified the effects of surface roughness, state changes (melt vs. solid), and particle size (e.g., Simurda et al. 2019;Williams and Ramsey 2019;Thompson and Ramsey 2020a). High temperature studies of cooling lava found that emissivity is significantly lower than previous assumed (Abtahi et al. 2002;Lee et al. 2013; Thompson and Ramsey 2020a,b). Therefore, it is likely that prior studies using a constant emissivity assumption overestimated heat flux and cooling rates and underestimated lava flow lengths (Harris et al. 1998;Harris and Rowland 2001;Harris 2013;Ramsey et al. 2019; Thompson and Ramsey 2020b).
Here, we establish the relationship between emissivity and temperature for basalt using ground-based multispectral TIR data acquired at a small lava flow emplaced in 2018 at Kīlauea volcano, Hawai'i. The temperature-dependent emissivity relationship we found became the basis for a new PyFLOWGO radiant heat flux module. The much larger channelized fissure 8 lava flow emplacement during the 2018 Lower East Rift Zone (LERZ) eruption was modeled using this new module to validate this approach. The modified model's effectiveness was assessed using high-resolution multispectral TIR data acquired on May 30, 2018.

Kīlauea volcano
Kīlauea is a basaltic shield volcano on the southeastern side of the Island of Hawai'i (HI, USA) that has erupted almost continually over the past 500 years. Mantle magmas are supplied to two reservoir systems below the main summit caldera at Kīlauea. The reservoirs are centered ~ 3-5 km beneath the south caldera and Keanakāko'i crater (Poland et al. 2014). From there, the magmas either ascend to a shallow reservoir and/or migrate throughout the > 80-km long rift system (Poland et al. 2014). Typically, recent eruptions have been effusive but phreatomagmatic explosive events have also occurred at the summit, which varied from short-lived to much longer (< 300 years) periods (Swanson et al. 2014). The effusive lava activity produces both 'a'ā and pāhoehoe (tubed and surface) flows with pāhoehoe flows commonly longer-lived Poland et al. 2014). Prior to May 2018, two main eruption styles were observed at Kīlauea volcano: (i) an overturning lava lake in the Halema'uma'u Crater from 2008 to 2018  and (ii) an extensive lava flow field emanating mainly from the Pu'u'Ō'ō vent (Wolfe et al. 1987;Heliker and Mattox 2003;Orr et al. 2013). More than 60 separate episodes of lava flows were observed over a 35-year period on the East Rift Zone with approximately 4.4 km 3 of lava emplaced, mostly as a series of tube-fed sheet-like and ropey pāhoehoe flows (Wolfe et al. 1987;Heliker and Mattox 2003;Orr et al. 2013;Neal et al. 2019).

Lower East Rift Zone eruption
In March 2018, the magma system beneath Kīlauea volcano started to pressurize, causing the lava levels at the Kīlauea summit and Pu'u'Ō'ō vent to increase (Neal et al. 2019). On April 30, 2018, ground deformation was observed down the rift system to the east and on May 3, 2018, lava erupted along a series of fissures near the Leilani Estates subdivision on the southeast of the island (Patrick et al. 2019). Over the next few weeks, lava erupted from a total of 24 fissures, with activity eventually focusing at the fissure 8 vent by May 27, 2018. Within 6 days, lava advanced 12.47 km to enter the Pacific Ocean at Kapoho Bay on June 3 ( Fig. 1a) (Neal et al. 2019). The eruption then remained centralized at fissure 8 for the next two months. From May 3 to May 27, the bulk effusion rates were estimated between 100 and 500 m 3 s −1 , but average rates increased to a maximum of ∼2000 m 3 s −1 reported in July (Neal et al. 2019;Patrick et al. 2019). Lava fountains ∼80 m high fed spatter that coalesced in a ∼30 m wide spillway before flowing into a ∼430 m wide perched channel ( Fig. 1b-c) (Patrick et al. 2019). Distally, the channel varied between ∼40 and ∼300 m wide. On August 4, effusion ceased after emplacing ∼1 km 3 of lava (Neal et al. 2019).

PyFLOWGO
FLOWGO is a one-dimensional numerical model that forecasts lava propagation confined within an open channel (Harris and Rowland 2001), which estimates the final flow length based on cooling. PyFLOWGO built upon the earlier FLOWGO model by improving initialization, iteration, and application criteria (Chevrel et al. 2018). The model combines the measured downflow changes in thermophysical properties and the best fit between results to produce a physically and thermally robust model (Harris and Rowland 2001). Downflow variations in velocity are calculated by estimating the crystallization, cooling, viscosity, and yield strength of the lava, with crystallization being the most influential factor. The model is coolinglimited with heat flux calculated until cooling increases the viscosity and reduces the velocity to a point where flow ceases (Fig. 2a) (Harris and Rowland 2001;Harris 2008). PyFLOWGO does not, however, replicate dispersive behavior observed at flow fronts because it models the flow in one dimension.
There are three main stages in the PyFLOWGO model: (1) determining the velocity of the lava based on a modified version of Jeffreys' equations; (2) calculating the heat flux from the lava during each propagation step to the next segment along the slope profile; and (3) determining the change in thermo-rheological conditions at each propagation step (Harris and Rowland 2001). Here, we address the third stage, which in turn impacts the other stages in the model.
Briefly, heat flux is originally calculated using a twocomponent temperature model of the lava surface representing the hot molten and cooler crust-covered end-members to determine an effective surface temperature ( T eff ) (Eq. 1) (Harris and Rowland 2001): where f crust is the fraction of surface crust cover, T c is the crust temperature, and T h is the temperature of the exposed molten flow. This effective temperature is then used to calculate heat flux (Eq. 3) and determine cooling rates (Harris and Rowland 2001). However, the original calculations do not account for the change in emissivity between these surfaces. A two-component emissivity model was developed to account for this change using an effective emissivity ( eff ) dependent on the fraction of the end-member surfaces (Eq. 2)  where c is the emissivity of the crust (0.95) and h is the emissivity of the hot molten surface (0.60) . The effective emissivity is combined with the effective temperature in the heat flux calculations (Eq. 3). This study improves on the two-component emissivity approach of Ramsey et al. (2019) by measuring the natural emissivity of crusted and molten surfaces and linking these to temperature and crust fraction. Our results show that it is important to represent the emissivity of both the crusted and molten lava surface because it has an influence on calculated cooling rates and ultimately, the modeled final flow length.

High-resolution ground-based data
The Miniature Multispectral Thermal Camera (MMT-Cam) acquires six, well-calibrated, surface radiance measurements between 7.5 and 12 µm at high spatial (< 1 m) and temporal (< 1 s) resolutions (Thompson et al. 2019). The data are calibrated for instrument attenuation and optical transmittance, as well as corrected for atmospheric emission and transmission effects using the SpectralCalc atmospheric simulator (Rothman et al. 2013;GATS 2019;Thompson et al. 2019). Surface radiance data are separated into surface kinetic temperature and emissivity using a modified Temperature-Emissivity Separation (TES) algorithm (Gillespie et al. 1998;Thompson et al. 2019). Total heat flux (ɸ tot ) (Eq. 3) and the fraction of exposed melt (Eq. 4) are then calculated where is the Stefan-Boltzmann constant, h c is the heat transfer coefficient, k is thermal conductivity, is thermal diffusivity, and A is pixel area. The heat transfer coefficient is calculated for a forced convection scenario because windy conditions were present during the acquisition periods (Harris 2013;Thompson and Ramsey 2020a). The fraction of exposed melt on the surface of the flow is calculated based on previous work deriving sub-pixel temperatures (Dozier 1981;Rothery et al. 1988;Harris 2013). In this study, we have adapted this to determine the faction of a pixel within a dataset that is completely molten and is similar to 1 − f crust in Eqs. (1) and (2), but derived directly from MMT-Cam data: where T liq is the liquidus temperature of a Hawaiian basalt (Abtahi et al. 2002).
On February 3, 2018, MMT-Cam data were acquired of a lava breakout from a tumulus on the coastal plain of Kīlauea volcano (Fig. 2b). We measured emissivity and temperature on the active region of the tumulus during the entire breakout event. The variability in both properties was analyzed and modeled for the six individual spectral bands and combined averages acquired by the MMT-Cam system. A sample was also collected to provide an estimate of vesicularity used in the model. Bulk density measurements revealed a bulk rock vesicularity of between 26 and 35% with a mean of 30%.
On May 30, 2018, the MMT-Cam was deployed on a helicopter over the LERZ to acquire data of the channelized flows originating from fissure 8 (Figs. 1, 2, and Suppl. 1). During the deployment, data were acquired of the lava fountain, spillway, perched channel, distal lava channels, and active flow front at the time (Fig. 1). Surface temperature, emissivity, fraction of exposed melt, heat flux, and channel widths were derived at each data acquisition location.

Topographic data
The pre-flow slope profile of the fissure 8 lava flow was extracted from the 10 m DEM derived from the United States Geological Survey's 7.5-min DEM Quads (National Oceanic and Atmospheric Administration 2007). The profile path was determined by overlaying a helicopter-based TIR map of Kilauea's LERZ fissure system produced on June 4, 2018 (U.S. Geological Survey 2018a). This was the day after the fissure 8 lava flow entered the Pacific Ocean at Kapoho Bay. The topography could have changed during the early phase of the eruption; however, the changes in slope profile were minimal. A syn-eruption DEM derived from an airborne Light Detection and Ranging (LiDAR) survey was conducted on July 8, 2018, until July 12, 2018 (U.S. Geological Survey 2018b). The dataset relieved a minimal (7.1% ± 4.1%) change in the channel slope profile as a result of the emplaced lava, therefore not significantly impacting the model interpretations.

PyFLOWGO modeling
The new temperature-dependent variable emissivity module developed from the MMT-Cam data was integrated into PyFLOWGO. The module was tested on the small tumulus-fed flow to constrain the new input parameters and analyze the sensitivity. The results were compared to the un-modified PyFLOWGO results using a constant emissivity value of 0.95 and 0.89, typical values used in previous thermal monitoring and modeling studies (e.g., Wright et al. 2008;Harris 2013;Patrick et al. 2017).
We first tested the updated PyFLOWGO model to simulate the 6.3-m tumulus-fed lava flow at 0.1-m intervals (model parameters in Table 1). The slope was calculated based on field measurements of the tumulus, with a change in elevation of 5.9 m over a horizontal distance of 10.0 m. Active channel widths were measured using the high spatial resolution (< 0.01 m) MMT-Cam data acquired during the emplacement. The initial channel width and depth measured at the bocca were 0.2 m and 0.3 m, respectively. These channel measurements, along with initial velocity measurements (~ 1.8 ms −1 ) acquired from the MMT-Cam data, were used to calculate an initial effusion rate of 0.11 m 3 s −1 . All other PyFLOWGO input parameters were either obtained from MMT-Cam data (e.g., crust cover fraction), through later analysis of samples (e.g., vesicularity), or based on assumptions using results from previous investigations on similar basaltic eruptions (e.g., viscosity) (Shaw 1972;Harris and Rowland 2001).
After testing, the modified model was used to simulate the 12.47-km-long lava flow that originated from fissure 8 at 10-m downflow intervals (model parameters in Table 2). The initial effusion rate (500 m 3 s −1 ) was constrained using reported values prior to the flow entering the Pacific Ocean on June 3, 2018 (Neal et al. 2019;Patrick et al. 2019). Channel width measurements were acquired periodically downflow using published TIR images (U.S. Geological Survey 2018a). On June 4, the average widths of the spillway, perched channel, and distal channel were 30 m, 340 m, and 150 m, respectively ( Fig. Suppl. 1). These values were used to validate the model results, but were acquired one day after the lava flow entered the Pacific Ocean, which may cause measurement discrepancies. PyFLOWGO calculates the channel width at each downflow propagation step. However, over time, channel widths can change, typically increasing after first forming. The initial viscosity was assumed to be 30 Pa·s to account for the increased lava mobility compared to the more viscous tumulus-fed lava flow (Shaw 1972;Harris and Rowland 2015;Dietterich et al. 2021). Table 1 The PyFLOWGO input models, modules, and parameters used to simulate the tumulus-fed lava flow 1 ∂φ/∂T cool = φ grown / (T erupt − T s ); 2 η melt = η erupt exp(0.04(T erupt − T core )); 3 η r = (1 − Rφ) −2.5 ; 4 τ 0 = ρ g h sin(α); 5 f crust = f init exp(α V mean ); 6calculated from MMT-Cam data at bocca using Eq. (4)

Models
Selection Reference

Tumulus-fed lava flow
The downflow evolution of two thermal properties (surface temperature and emissivity) and two thermophysical parameters (fraction of exposed melt and heat flux) was observed during the tumulus-fed lava emplacement over a 40-min period (Fig. 3). The MMT-Cam data transects show high initial temperatures (up to 1450 K) 1-2 m and 3-4 m from the bocca, where the largest active breakouts were observed. As the emplacement progressed, the highest temperatures progressed farther from the bocca, with major peaks at 3-4 m, ∼4.7 m, and 5.2-6.3 m (Fig. 3a). Lower emissivity values (> 0.65) correlated with these highest temperatures, especially at 3.5-3.8 m and 5.5-6.3 m (Fig. 3b). The fraction of exposed melt is similar  . 3 Temporal transects of the a surface temperature, b six-point average emissivity, c fraction of exposed melt, and d heat flux along the central channel of the 6.3-m tumulus-fed lava flow. The original bocca is at 0 m (see Fig. 2b). The MMT-Cam data were acquired at an approximate distance of 10 m from the target. The transect line color transition from purple to yellow with time. The spatial resolution is ~ 0.015 m to the temperature with an aerial exposed melt fraction near 0.99 between 1-2 and 3-4 m and a minima of 0.15 in the less active regions. Higher fractions of exposed melt were observed farther from the bocca as time progressed. Heat flux data were similarly variable, with maximum values of 0.18 MW observed (Fig. 3d). Overall, the highest thermal values (lowest emissivity) transition farther from the bocca with time (each transect line in Fig. 3 from purple to yellow, at ~ 0.5-s intervals).

fissure 8 lava flow
The helicopter-borne MMT-Cam acquired data at three locations along the fissure 8 lava channel. The targeted sites were the ∼80-m-high lava fountain at the vent, the perched lava channel, and the distal lava channels. At the vent, the highest temperatures, fraction of exposed melt, and total heat flux were observed at the base of the fountain, where the emissivity values are the lowest (Fig. 4). The fraction of exposed melt clearly shows the pathway of molten lava. The TIR data also highlighted smaller breakouts of molten lava to the south of the vent, away from the main channel. Generally, the average emissivity was the lowest around the vent and in some channel pathways in the spillway. It was also lower above the lava fountain likely due to the absorption of SO 2 .
MMT-Cam data of the perched lava channel revealed a complex surface of molten lava and rafted crust (Fig. 5). Rafted plates dominated the surface of the channel (> 80%). The average emissivity was moderate, with the majority of the channel surface having a value of ∼0.75. Ribbon-like channels of lower emissivity (as low as 0.6) correlate with areas of molten lava. Heat flux values were lower in the perched channel compared to the vent region; however, elevated flux was detected in small areas suggesting the presence of subsurface lava pathways. Overall, the thermal properties of the perched lava channel indicate a fairly well-insulated flow. The distal lava channel was comprised of diverse morphologies with variable thermal properties, including poorly to well-insulated rafted portions, and small roofed over segments (Fig. 6). The surface temperature and fraction of exposed melt data distinguish between cooler crusted lava rafts and exposed molten regions. Surface temperatures up to 1475 K were measured on the molten surfaces, compared to average temperature of ∼1000 K on the crusted rafts. A similar difference between molten and crusted regions was observed in the fraction of exposed melt, with average calculated values of ∼0.90 and ∼0.55, respectively. The cooler crusted lava rafts were constrained to the center of the channels whereas the molten regions were located along channel margins. Not surprisingly, the greatest heat flux was observed in the main channel along those margins. The lateral variability in thermal properties was the greatest in wider channels. Margins presented the lowest emissivity regions, with values as low as ∼0.6, suggesting the presence of a complex surface of cooling lava and minor breakout/overflow events.

Variable emissivity module
MMT-Cam measurements of the tumulus-fed lava were first performed to establish the relationship between emissivity and temperature for a basaltic melt (~ 50 wt.% SiO 2 ). Lava temperatures ranged from 900 to 1450 K (Figs. 7a and Suppl. 2). The 8. 04, 8.55, 9.55, and 10.04 µm bands revealed a strong inverse correlation between surface temperature and emissivity. There was a minor inverse relationship in the 8.99 µm band and a positive relationship in the 11.35 µm band (Suppl. Fig. 2). Additionally, the absorption shallows and broadens during cooling (Fig. 7a). Trends were then combined with equal weighting (as the bands have the same full width half maximum) to characterize the average  (Fig. 7b).
The emissivity relationships established with MMT-Cam data were used to modify the radiant heat flux calculations of PyFLOWGO, specifically to vary emissivity with lava surface temperature. Regression analysis indicated that both linear and quadratic models accurately represent the relationship between surface temperature and emissivity of lavas during cooling (Figs. 7b and Suppl. 2 and Suppl. Table 1). The linear regression was chosen for module development into the PyFLOWGO model because of its relative simplicity and higher coefficient of determination (∼0.94). The average 6-point emissivity values were used, as no perceivable difference was achieved using the more complex, separate emissivity band regressions.

Tumulus-fed lava flow
PyFLOWGO modeling revealed that the final f low length increased by ∼5% using the new variable emissivity module compared to the previous constant emissivity (0.95) module. Importantly, the longer f low length better matched the field observations (Figs. 8 and 9). The increased length is a result of lower heat flux (reduced cooling rate), which maintains higher temperature and therefore lower viscosities and crystallization, and, ultimately, delays stopping (Figs. 8b, d, e, and f). The greatest impact was observed within the final 20% of the flow. For example, at 6.0 m, where the constant emissivity (0.95) module simulated the flow stalling, the variable emissivity module estimated ∼75% less heat flux and 4.9-5.4% higher temperatures (Figs. 8b, d-f).
The new variable emissivity module for PyFLOWGO was validated by comparing the modeled channel widths, surface temperature, and total heat flux results with those derived from MMT-Cam data (Figs. 8c-e). The modeled channel widths show a good correlation with MMT-Cam results until the flow reached a length of ∼4.5-5 m. Beyond 4.5 m, the tumulus-fed lava flow separated into multiple smaller channels causing the width of the main channel to decrease, something that is not simulated by PyFLOWGO. The field surface temperature measurements show a wide variation due to the small scale (< 0.1 m) changes in surface cooling downflow that is likely missed in models. The total heat flux calculated in both the updated and unmodified PyFLOWGO model was < 100% higher than those derived from the MMT-Cam data (Fig. 8e). The disagreement in values increased with flow length and was over an order of magnitude different farther downflow.
In these PyFLOWGO model simulations, there were minimal differences between the emissivity modules for the initial ∼50% of the flow length (Fig. 9). The greatest difference was observed in the final ∼20% with the variable emissivity module producing a longer flow. Viscosity validation calculations are within an order of magnitude of the model results. Generally, the rates of change in values of thermo-rheological properties (e.g., Δv/Δx) using this variable emissivity module were comparatively slower with propagation and flow evolution. For example, the rate of change in mean velocity is ∼3.2% slower, which provides improved insights into when a lava flow may reach a downflow location (Fig. 9d). Similar trends were observed in the other properties (Figs. 9a-c).

fissure 8 lava flow emplacement
The variable emissivity module resulted in a better match to the actual length of the fissure 8 lava flow, increasing by ∼7%, from 10.68 to 12.47 km (Figs. 10 and 11). Similar to the tumulus modeling, the increase was caused by the reduction in cooling rates as a result of lower calculated radiant heat flux (∼32%). The variable emissivity module also simulated comparatively higher temperatures (∼9%) throughout the flow length, as a result of slower cooling rates. The differences between the two emissivity modules increased progressively with distance (Fig. 10).
The variable emissivity PyFLOWGO results of the fissure 8 lava flow emplacement were validated by comparing the modeled channel widths, surface temperatures, and heat flux with field measurements (Fig. 10 and Suppl. 1). There is strong agreement between the model results and the field measurements, with the variable emissivity module having a better agreement at greater flow lengths (> 9 km). The greatest difference was observed in the proximal region (0.1-4.0 km from the vent) where the wide perched lava channel formed. The thermo-rheological results calculated using the variable emissivity module were similar to the original emissivity module but were more consistent with slower cooling rates associated with less efficient radiant heat flux (Fig. 11). This caused lower crust cover fractions and viscosities, as well as higher core temperatures and mean velocities to be simulated with distance, compared to  Fig. 2b). The gray region highlights the 2-sigma variance of emissivity and the black line represents the mean. The orange and green lines represent the computed linear and quadratic regressions, respectively. Note, only a random 50% of the data are plotted here due to graphic limitations, but all the data are used in the calculations.
previous modeling efforts. Overall, the average rate change of these thermo-rheological properties was ∼7% lower using the variable emissivity module and the ocean entry length was more consistent with the field measurements for all the results (12.47 km).

Discussion
This study found that the modeled heat flux is overestimated compared to field measurements (Fig. 8), similar to previous studies ). This could be attributed to two possibilities. First, the across-flow thermal variability is poorly constrained by the model. Typically, there is a nonlinear temperature gradient resulting in non-uniform heat flux across the flow channel. The discretized calculation in PyFLOWGO oversimplifies this variability by assigning a singular thermal component to the entire width for each propagation increment (e.g., every 0.1 m for the tumulusfed flow) based on the linear two-component thermal mixing calculation (Eq. 1). The MMT-Cam has a higher spatial resolution (< 0.01 m) both across and downflow, capturing the spatial complexity of the thermal variability. The heat flux is calculated for each image pixel and all the values summed over a comparable area to the model. Furthermore, as the channel width increases, the model-calculated heat Fig. 8 Comparison between the a emissivity, b effective temperature, c channel width, d surface temperature, e total heat flux, and f lava core temperature produced by the PyFLOWGO model during the tumulus-fed lava flow. These were computed using the variable emissivity module (black solid line) and constant emissivity module (blue and orange dashed line). Input parameters are in Table 1. Red vertical line represents the measured final flow length. Field measurements (red dots) derived from MMT-Cam data are used to validate channel widths, surface temperatures, and total heat flux results flux will become increasingly overestimated. Second, the air directly above the flow is heated producing a thermal gradient between the camera and the surface. This gradient, if not accurately corrected, may lead to an inaccurate heat flux attributed to the lava surface.
Channel widths are incorrectly modeled in two locations. PyFLOWGO overestimates the channel width in the final ∼20% of the flow due to the loss of a main channel and divergence of the lava into multiple smaller channels. The model becomes poorly constrained in this situation (Harris and Rowland 2001). The agreement was improved by combining the widths of the multiple channels observed at the flow fronts in the TIR data. This increased the final channel widths for the tumulus-fed and fissure 8 lava flow emplacements to ~ 2 m and ~ 1800 m, respectively. The fissure 8 lava flow emplacement simulations also underestimate the width of the proximal, well-insulated, perched lava channel because the model is not designed to simulate well-insulated roofed channels (see Patrick et al. 2019;Figs. 1b, 5, and 10b). Nevertheless, PyFLOWGO was able to accurately model the exposed channel emanating from this region and the ocean entry length, implying a completely exposed channelized flow is not required for accurate modeling.
A similar relationship is observed between the model results and field measurements of surface temperature and viscosity (Figs. 8, 9 and 10). The field measurements vary randomly compared to the model results. This variability is a result of the complex nature of lava surfaces during propagation and cooling, the small-scale interaction with pre-existing topography, and the discrete sampling bias of the field measurements. Nevertheless, the field measurements do constrain these results in the model with the surface temperature Fig. 9 Thermo-rheological variations of the tumulus-fed lava flow comparing the PyFLOWGO model simulation results using the variable emissivity (black solid line) and constant emissivity (blue and orange dashed line) modules. The results using the original emissivity module consistently under predicted the actual flow length (red vertical line). Field measurements (red dots) were derived from MMT-Cam and geometry data of the tumulus-fed lava flow of viscosity using a modified Jeffreys' equation (Nichols 1939) and viscosity results within 25% and an order of magnitude, respectively. The fissure 8 lava flow modeled velocity results near the vent and along the proximal lava spillway agree with other velocity measurements taken in the field (Patrick et al. 2019;Dietterich et al. 2021). Proximally (< 2 km), the modeled velocities oscillated between ~ 10 and 28 ms −1 with other field studies measuring velocity oscillating between ~ 5 and 20 ms −1 , providing validation for the fissure 8 lava flow emplacement model results.
During the development of the variable emissivity module, scatter in the MMT-Cam temperature and emissivity data led to some uncertainty in the calculated linear regression (Fig. 7b). The effect of this uncertainty on the PyFLOWGO results was determined using a Monte Carlolike methodology. This constrained the variability in the final flow length by randomly repeating the flow propagation simulations 10,000 times (Fig. 12). The regression constants were randomly generated over the 2-sigma standard deviation data range in a normal Gaussian distribution (Figs. 12b and d). The uncertainty in the final length of the tumulus-fed lava flow was ∼0.5% with a total variability of ∼3% (Fig. 12a), whereas the uncertainty in the final length of the fissure 8 lava flow was ∼2% with a total variability of < 1% (Fig. 12c). Therefore, uncertainty due to the linear The greatest deviation in emissivity between the modules was observed between the lava liquidus (∼1425 K) and solidus (900-1000 K) temperatures (Cashman et al. 1999;Gottsmann et al. 1999), causing a similar trend in the thermo-rheological properties (Fig. 11). Viscosity has a strong control on the other properties in the model and flow length as it is calculated based on the derived core temperature using Dragoni's (1989) equation. The greatest difference in derived properties calculated by using the variable emissivity module in PyFLOWGO was observed farther downflow from the vent/bocca. The difference between the two modules (variable and constant emissivity of 0.95) became most significant (> 30%) where the core temperature decreased below 1235 K (e.g., at ∼10.5 km in the fissure 8 lava emplacement) (Figs. 9 and 11). The greatest increase in emissivity is observed at lava temperatures lower than the liquidus as a result of crystallization and surface crust formation, including viscoelastic and glassy crusts (Fig. 7) (Thompson and Ramsey 2020b). This is diagnostically observed in the emissivity data by a broadening and shallowing of the absorption features (e.g., Fig. 7a), as a result of molecular structural and vibrational changes during cooling and crystal/crust formation (Grzechnik and McMillan 1998;Lee et al. 2013).
As a result of this study and implementation of the variable emissivity module in PyFLOWGO, cooling rates are better constrained through improved radiant heat flux, crust formation, and viscosity calculations. The variable emissivity relationship for calculating radiant heat flux within PyFLOWGO improves the constraint on other thermo-rheological/physical properties derived by the model. The variable emissivity relationships more accurately represent the thermal changes observed in nature over the use of single emissivity value, resulting in a ~ 2 km increase in the flow length prior to the ocean entry (for a flow ~ 12 km long). If the variable emissivity module was applied to longer subaerial lava flows at Etna or Piton de la Fournaise, models using this variable emissivity model could predict final flow lengths several kilometers longer. These extra kilometers could cause large socioeconomic impacts to these areas, including to residences, industries, and transportation networks. Therefore, it is critical to accurately model lava flows in these regions to provide hazard response agencies with as much reliable information as possible to ultimately reduce the risk lava flows pose to local societies.
However, this new temperature-dependent emissivity module is only applicable for lava with similar eruption and liquidus temperatures as the Hawaiian basalts modeled in this investigation. For example, lavas with different silica and alkali contents will have unique eruption and liquidus temperatures. Therefore, these would have a different temperature-dependent emissivity variability (Lee et al. 2013). As a result, similar investigations need to be conducted to evaluate the variability in emissivity with temperature for a variety of lava compositions. Additionally, the initial viscosity estimate has a limited influence on the final flow length. For example, if the initial viscosity was increased from 30 to 200 Pa s in the fissure 8 lava flow emplacement model, the ocean entry PyFLOWGO 10,000 times. The emissivity distributions (a and c) represent the average emissivity computed during the entire modeling duration. In comparison, the flow lengths achieved using the constant emissivity of 0.95 were b 6.0 m and d 10.68 km length would increase by ~ 100 m (~ 0.7%). This is compared to the ~ 2 km increase in final flow length achieved by using the variable emissivity module. Nevertheless, more studies are required to evaluate these and other thermophysical parameters influence on models. These investigations could be accomplished in a laboratory or the field to enable similar lava flow propagation modeling, developed in this study, to be applied to more lava flows in near real-time.
This study has highlighted multiple potential avenues for future investigation in lava flow propagation modeling and the data required to improve these models. Investigations are required to increase the ability to model the across-flow variability and accurately adapt these measurements during propagation. Additionally, we need to increase our understanding of the petrological implications of varying emissivity in these models to improve the physical and rheological calculations. An increase in high-spatiotemporal (< 10 mm spatial and < 1-s temporal resolution) multispectral (> 5 spectral bands) TIR data is required to conduct these investigations in situ and ultimately further improve lava flow modeling to reduce the risks posed to nearby populations.

Conclusion
Emissivity is an important material property that must be considered in calculations of radiative cooling experienced by molten lava during propagation and cooling. The efficiency by which a surface cools radiatively is strongly controlled by the physical state of the surface and consequently the surface temperature, as well as other properties (e.g., composition and surface morphology). Variable emissivity is shown to have a measureable effect on heat flux (> 30% decrease), which translates to the potential of longer flows (~ 7% increase). Lava propagation modeling assuming a constant emissivity value close to 1.0 in heat flux calculations will underestimate the final flow length resulting in inaccurate lava flow hazard potential.
Here, we used high-resolution multispectral TIR data of a tumulus-fed lava flow on Kīlauea volcano to derive the complete thermal evolution and the correlation between emissivity and temperature during propagation. The relationship was used to develop a new variable emissivity module that was integrated into the PyFLOWGO model. This improved the accuracy of calculated thermo-rheological properties during flow propagation and cooling. The modified PyFLOWGO model was used to simulate the 2018 fissure 8 lava flow emplacement at Kīlauea volcano. These results were compared with those derived using the original PyFLOWGO model and validated using field observations. The variable emissivity module produced a lower final heat flux (~ 75%) and a more accurate length fit (increased by ~ 7%; ~ 2 km) compared to using a constant emissivity.
This study has shown that PyFLOWGO already models lava propagation quite well when the input parameters are well known and well constrained. As such, the model is highly adaptive, having been applied to a variety of flows on Earth and other terrestrial bodies in the past (e.g., Rowland et al. 2004;Harris et al. 2019;Ramsey et al. 2019). All lava flow modeling is ultimately limited by the accuracy of input parameters and prior thermo-rheological measurements of the lava, as well as the ability to accurately simulate the physical and chemical processes occurring as the lava propagates. The modified PyFLOWGO model developed here, using a new variable emissivity module developed from unique in situ measurements of basalt temperature and emissivity, increased the model's accuracy. It is hoped that the development of this new module will increase its accuracy and reduce the processing time in future eruptions, which will better inform hazard assessments and reduce the vulnerability of local populations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.