Linking lava morphologies to effusion rates for the 2014–2015 Holuhraun lava flow field, Iceland

Determining the parameters that control fissure-fed lava morphologies is critical for re-constructing the complex emplacement histories of eruptions on Earth and other planetary bodies. We used a geomorphological map of the 2014–2015 Holuhraun lava flow field, in combination with new constraints on lava emplacement chronology and two independently derived time-averaged discharge rate (TADR) data sets, to analyze correlations between lava morphology and effusion rate. Results show that lava morphologies are dominantly controlled by effusion rate at the vent during the early phases of the eruption and by lava transport processes as the system evolves. Initially, TADR and its variance, which reflect pulsation in the lava supply rate from the vent, directly affect lava emplacement styles. However, as the eruption progresses, the lava transport system exerts a stronger control with channels and ponds that can either dampen variation in local effusion rate or create surges during sudden drainage events. The Holuhraun eruption predominantly produced rubbly lava in its earlier eruption phases and transitioned into the production of spiny lava toward the end of the eruption. However, a drop of TADR during the first phase of the eruption correlates with a decrease in rubbly lava formation and an increase in spiny lava production. This suggests that a change in effusion rate caused the observed transition in lava type. Our findings show that rubbly lava is formed under relatively high local effusion rates with pulsating supply conditions, whereas spiny lava is formed under lower local effusion rates and steadier supply.


INTRODUCTION
Linking lava morphologies to physical parameters is key for understanding how post emplacement characteristics can be used to reconstruct eruption conditions on Earth and other planetary bodies-especially where lava flow fields record a complex history of initial flow emplacement followed by processes of dis ruption and inflation.
Many parameters have been suggested as controls for lava morphologies (see the Supplemental Material 1 ), but, generally, mor phologies are affected by a combination of the initial chemical and thermophysical properties of the lava (Macdonald, 1953;Williams and McBirney, 1979) and dynamic processes, such as volumetric flow rate and shear strain rate (Rowland and Walker, 1990;JuradoChichay and Rowland, 1995), which in turn depend pri marily on flow thickness and viscosity, local slope, gravitational acceleration, and the exter nal environment (Glaze et al., 2014). However, interconnections between these parameters are complex, and while previous studies (e.g., Peter son and Tilling, 1980;Hon et al., 2003;Chevrel et al., 2018) have greatly enhanced our under standing of the relationships among lava shear strain rate, viscosity, and the development of pāhoehoe and ʻaʻā, the factors controlling the development of transitional lava types remain poorly constrained.
Here, we test the hypothesis that for basaltic fissurefed eruptions of constant composition, effusion rate and its variance are the principal controls on the different lava morphologies. To accomplish this goal, we created a novel chronological map of the flow growth for the 2014-2015 Holuhraun lava flow field in Iceland, and investigated relationships between timeav erage discharge rates (TADRs; Coppola et al., 2017;Bonny et al., 2018) and observed lava morphologies, or facies (Voigt et al., 2021). The 2014-2015 Holuhraun event provides an excel lent case study because it was well monitored with remotesensing data and other timeseries observations (Gudmundsson et al., 2016;Peder sen et al., 2017), its lava composition remained essentially constant (Halldórsson et al., 2018), the lava was emplaced on a nearly uniform gla cial outwash plain with an average slope of less than 0.3° (Voigt et al., 2021), and most flow units remained well exposed after the end of the (6 month) eruption.

HOLUHRAUN GROWTH, FACIES, AND TIME-AVERAGED DISCHARGE RATE
The 2014-2015 Holuhraun eruption ex truded a denserock equivalent volume of ∼1.21 km 3 (Bonny et al., 2018) and covered 83.82 km 2 (Voigt et al., 2021) between 29 Au gust 2014 and 27 February 2015. To map the growth of the lava field (Fig. 1A), we combined GPS measurements of lava margins acquired in the field (Pedersen et al., 2017) with airborne radar data (Pedersen et al., 2017) and satellite radar data from TerraSARX and TanDEMX (Dirscherl and Rossi, 2018), LANDSAT7 and LANDSAT8, EO1 Advanced Land Imager (ALI) and Hyperion, and SENTINEL1, and Moderate Resolution Imaging Spectroradiom eter (MODIS) and National Oceanic and Atmo spheric Administration (NOAA) Advanced Very High Resolution Radiometer (AVHRR) images (Jónsdóttir et al., 2014). Our chronological map contained 49 units and thus included enough time steps to test our hypothesis.
Holuhraun lava facies are defined by re gions of similar albedo and texture mapped at 1:800 scale, with unit names corresponding to the dominant lava type within each unit (Voigt et al., 2021). We focused on the three most extensive facies, which together comprise 93% of the lava flow field, namely, the rubbly (57%), spiny (26%), and undifferentiated rubbly-spiny (10%) facies. From the facies and chronologi cal maps, we calculated the area of each facies emplaced per chronological unit and normal ized that by the total area of all facies emplaced within the same chronological unit (Fig. 2). The area percentage emplaced by a facies per day was then retrieved using linear interpolation (see the Supplemental Material).
TADR is defined as the volume flux aver aged over a specific time, whereas local effusion rate is defined as the instantaneous effusion rate for a specific flow unit (Harris et al., 2007). To date, two different TADR data sets are available for the Holuhraun eruption: (1) satellitederived TADR compiled by Coppola et al. (2017); and (2) a combination of groundbased and satel litederived TADRs compiled by Bonny et al. (2018). In this work, we refer to values provided by Coppola et al. (2017) in the text (unless oth erwise noted), but we show results of both input data sets in the Supplemental Material. The two TADRs follow the same general trend, but the magnitudes of the values differ. Because we are primarily interested in analyzing the correlation, the overall TADR patterns are more important than their absolute values.

CONTROL OF LAVA FACIES BY THE EFFUSION RATE Hypothesis
We hypothesize that lava emplacement conditions depend upon local effusion rate, but assume that these local conditions inherit a damped, but detectable signal of source vent activity. We further assume that the TADR pro vides an adequate proxy for the effusion rate at the vent and that the TADR primarily reflects variations in lava supply conditions. However, we acknowledge that the lava transport sys tem modifies local effusion rates and that flow branching and sudden release of lava stored within reservoirs, respectively, would either decrease or increase the local supply of lava. Under these assumptions, our central hypoth esis is that the formation of rubbly lava is fa vored under pulsating effusion rate conditions (reflected in high TADR variance) because the initially coherent crust would be fragmented into pieces. Additionally, we expect that the forma tion of rubbly lava would be favored when the local effusion rate is high, leading to a positive correlation between the TADR and the emplace ment of rubbly facies.
In contrast, the spiny facies is hypothesized to form under steadier lava supply conditions, reflected in a lower local effusion rate with smaller variance. This is because the spiny fa cies typically exhibits a coherent crust. Spiny facies emplacement is therefore expected to be anticorrelated with the TADR and its variance.
The undifferentiated rubbly-spiny facies is a hybrid between the two principal facies: rubbly and spiny. This facies exhibits highly fragmented surfaces, but it also has regions of coherent crust, which are exposed in close spa tial association to one another. We expected this facies to form when conditions favorable for the formation of rubbly lava changed to conditions favoring spiny lava, or vice versa.

Statistical Evidence
In general, both TADR data sets (Coppola et al., 2017;Bonny et al., 2018) show a decreas ing trend during the course of the eruption, and therefore we expect that there would be a change in the predominant facies emplaced per day over time. Our results show that there were three phases in which one facies dominated (Fig. 2). During phase I (this work defines day 0 to be 29 August 2014 at 12 p.m. GMT), which lasted until day 77, the rubbly facies was the dominant eruption product. In phase II, from day ∼77 to 91, undifferentiated rubbly-spiny facies domi nated. In phase III, starting at day 91, the lava flow field primarily grew by the emplacement of spiny facies. Figure 2 also shows a temporary drop in TADR in phase I starting on day 17 in Bonny et al. (2018) data and day 19 in Coppola et al. (2017) data. During this time, the dominant facies changed from the rubbly facies at high TADRs (peaks up to ∼419 m 3 /s on day 13) to spiny facies at relatively low TADRs (78 m 3 /s on day 19). When the TADR increased, around day 24, the dominant facies switched back to rubbly (∼250 m 3 /s, and even 440 m 3 /s in Bonny et al. [2018] data). Between days 25 and 28, the TADR decreased, and spiny forma tion dominated, and then gradually increased to 80%-85%. This shows that lava facies correlate with the TADR and is not simply a monotonic change from one facies to another with time.
The initial period of rubbly facies emplace ment closely followed the TADR until day ∼33 (Bonny et al., 2018), or ∼47 (Coppola et al., 2017). Between days 25 and 80, the peaks of the TADR correlate well with the peaks in rubbly facies production; however, the absolute values of the TADR are smaller compared to the first month, ranging from 120 to 320 m 3 /s in Coppola et al.'s (2017)  In the earlier stages of the eruption, the TADR variance was highest until day 35 in the Bonny et al. (2018) data. In the Coppola et al. (2017) data, this period variance extended up to day 45. Thereafter, the TADR decreased steadily in both data sets. During phase II, the TADR variance was close to zero in the Bonny et al. (2018) data and varied only minimally in the Coppola et al. (2017) data until day 91. In phase III, the vari ance was close to zero in both data sets.
To quantify the correlations, we calculated the Spearman correlation coefficient, r s (Spear man, 1904). We chose the Spearman coefficient because the correlation of our data was expected to be monotonic but not linear. In contrast, the Pearson correlation coefficient assumes a linear relationship of the quantities to be correlated. Table 1 summarizes the resulting r s values for the correlation of TADR and facies emplaced per day.
The rubbly facies showed a strong correla tion of 71% with the TADR normalized per day and a correlation of 61% with the variance. In contrast, the spiny facies showed an anticorre lation of -66% with the TADR, as well as an anticorrelation of -61% with its variance. The undifferentiated rubbly-spiny facies had moder ate correlations of 54% with the TADR and of 58% with its variance.
On day 36, growth stopped on the north ern part of the lava flow field and a new flow unit developed focused in the southern portion (Fig. 1). This lava unit was confined by a former river channel. Around day 44, a new lava unit developed, hosting a lava pond. This lava pond acted as a major distributor for the flow units emplaced until day 94. The pond stored lava and thus acted as a buffer, resulting in a rapidly decreasing TADR variance. Around day 60, the pond size started to decrease, but it was still feeding rubbly lobes to the south, resulting in 75%-100% rubbly formation until day 77.
At the end of phase I (Fig. 2), the emplace ment of rubbly lava stopped until a second period of rubbly lava formation began on day 105, when new lava lobes developed northwest of Baugur (i.e., the main vent for the eruption; Figs. 1 and 3), starting with a rubbly production rate of 40% per day and subsequently increas ing to ∼95% and 80% on eruption days 129 and 131, respectively. The rubbly emplacement during this time was not correlated with a high TADR, nor with a high variance of TADR.
These results suggest that in the early stages of the eruption, the lava output at the vent itself (i.e., TADR for the eruption at a specific time) was equivalent to the local effusion rate, and we obtained a very high correlation between TADR and emplaced facies. In the case of variations in the local effusion rate, due to changes in the lava transport system, such as the development of the lava pond, which can cause a sudden surge in lava flux within a newly formed pathway, the  A C B TADR is not a good proxy for the local effu sion rate. Instead, these cases require additional constraints from other remotesensing and field based observations to constrain the timing and circumstances of rubbly lava breakouts. Evidence of the effects of sudden changes in local effusion rate is available for the lobe that began its emplacement on day 111 onto the gla cial outwash plain (sandur) northwest of Baugur ( Figs. 1 and 3). The northern levée of one of the main channels (Fig. 3) collapsed, spawning a breakout to the north, which formed an active lava channel. This event implies that the lava lobes were fed with high local effusion rates, which is consistent with surges in rubbly facies growth on days 129 and 131. On day 135, the lava channel that was feeding the rubbly lobes in the northwest of the lava flow field still carried incandescent lava, but then started to roof over. On day 145, the entire channel developed a crust and was cooling with no exposures of incandes cent lava on the surface. The flow field growth with rubbly lava ended after the emplacement of the northwestern lobes around eruption day 134.

CONCLUSION
Analysis of the lava morphology (i.e., facies) combined with the chronological emplacement of lava units and TADRs for the 2014-2015 Holuhraun eruption in Iceland demonstrates a strong correlation between TADR and the final lava products. However, our findings indicate that it is most likely not just the TADR, but rather the local effusion rate that exerts the dominant control on the morphologies of transitional lava types. Our results suggest that rubbly facies forms at high local effusion rates with stronger pulsating lava supply, and spiny facies forms un der a low local effusion rate with a rather con tinuous lava effusion. Generally, the two TADR data sets correlate well with the final lava mor phology (ranging from a strong correlation of 71% to moderate correlations of 54%). In the early stages of the eruption, when TADR and effusion rate are similar by definition, the corre lation is even stronger. In the later stages, rubbly facies was emplaced under neither high TADR nor high variance conditions; however, qualita tive evidence for high local effusion rates exists, explained by a levée breaching event, causing a temporary surge in the local effusion rate. This information is particularly important for interpreting eruption conditions for ancient lava flow fields on Earth and other planetary bodies where only a postemplacement geologic record is available. Regions like Elysium Planitia on Mars show lava morphologies with very similar characteristics to the Holuhraun eruption prod ucts, with both highly disrupted and continuous surfaces that exhibit inflation features (e.g., Kes zthelyi et al., 2004;Vaucher et al., 2009;Jaeger et al., 2010;Voigt and Hamilton, 2018). The spa tial relationship of these two facies is very similar to that of the rubbly and spiny facies within the Holuhraun lava flow field. This implies that both Holuhraun and analogous lava flows on Mars were dominated by an initial lava emplacement phase with relatively high local effusion rate and variance, which then transitioned into the em placement of spiny lava as local effusion rates decreased to become steadier and more continu ous. Thus, critical information about effusion rate can be revealed for past fissurefed eruptions by carefully examining the lava morphologies preserved within the geologic record.

ACKNOWLEDGMENTS
We thank Oryaëlle Chevrel, Scott Rowland, and one anonymous reviewer for their thoughtful reviews, which greatly improved this manuscript. J.R.C. Voigt thanks the Geological Society of America (GSA) for a 2016 GSA Lipman Research Award, and the NASA Future Investigators in NASA Earth and Space Science and Technology program grant NNH19ZDA001N FINESST; C.W. Hamilton acknowledges NASA support from Planetary Science and Technology from Analog Research grant 80NSSC21K0011 and the Fulbright-National Science Foundation Arctic Research Scholarship program, administered by Ful bright-Iceland. We thank Diego Coppola for providing the timeaveraged discharge rate (TADR) data, and Ulrich Münzer and Iceland subglacial volcanoes inter disciplinary early warning system (IsViews) for the TerraSARX, TanDEMX data, and UltraCam data.