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

Linking lava morphologies to physical parameters is key for understanding how postemplacement 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 disruption and inflation.

Many parameters have been suggested as controls for lava morphologies (see the Supplemental Material¹), but, generally, morphologies 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; Jurado-Chichay and Rowland, 1995), which in turn depend primarily on flow thickness and viscosity, local slope, gravitational acceleration, and the external environment (Glaze et al., 2014). However, interconnections between these parameters are complex, and while previous studies (e.g., Peterson and Tilling, 1980; Hon et al., 2003; Chevrel et al., 2018) have greatly enhanced our understanding 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 fissure-fed 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 time-average 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 excellent case study because it was well monitored with remote-sensing data and other time-series observations (Gudmundsson et al., 2016; Pedersen et al., 2017), its lava composition remained essentially constant (Halldórsson et al., 2018), the lava was emplaced on a nearly uniform glacial 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.

The 2014–2015 Holuhraun eruption extruded a dense-rock equivalent volume of ~1.21 km³ (Bonny et al., 2018) and covered 83.82 km² (Voigt et al., 2021) between 29 August 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 TerraSAR-X and TanDEM-X (Dirscherl and Rossi, 2018), LANDSAT-7 and LANDSAT-8, EO-1 Advanced Land Imager (ALI) and Hyperion, and SENTINEL-1, and Moderate Resolution Imaging Spectroradiometer (MODIS) and National Oceanic and Atmospheric 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 regions 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 chronological maps, we calculated the area of each facies emplaced per chronological unit and normalized 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 averaged 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) satellite-derived TADR compiled by Coppola et al. (2017); and (2) a combination of ground-based and satellite-derived TADRs compiled by Bonny et al. (2018). In this work, we refer to values provided by Coppola et al. (2017) in the text (unless otherwise 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.

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 provides 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 system 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 hypothesis is that the formation of rubbly lava is favored 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 formation of rubbly lava would be favored when the local effusion rate is high, leading to a positive correlation between the TADR and the emplacement 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 facies 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 spatial 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.

In general, both TADR data sets (Coppola et al., 2017; Bonny et al., 2018) show a decreasing 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 dominated. 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³/s on day 13) to spiny facies at relatively low TADRs (78 m³/s on day 19). When the TADR increased, around day 24, the dominant facies switched back to rubbly (~250 m³/s, and even 440 m³/s in Bonny et al. [2018] data). Between days 25 and 28, the TADR decreased, and spiny formation 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 emplacement 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³/s in Coppola et al.'s (2017) data, and even lower in Bonny et al.'s (2018) data.

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 variance was close to zero in both data sets.

To quantify the correlations, we calculated the Spearman correlation coefficient, r_s (Spearman, 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 correlation of 71% with the TADR normalized per day and a correlation of 61% with the variance. In contrast, the spiny facies showed an anticorrelation of–66% with the TADR, as well as an anticorrelation of–61% with its variance. The undifferentiated rubbly–spiny facies had moderate correlations of 54% with the TADR and of 58% with its variance.

On day 36, growth stopped on the northern 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 emplacement 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 increasing 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 TADR is not a good proxy for the local effusion rate. Instead, these cases require additional constraints from other remote-sensing 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 glacial 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 incandescent lava on the surface. The flow field growth with rubbly lava ended after the emplacement of the northwestern lobes around eruption day 134.

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 under a low local effusion rate with a rather continuous lava effusion. Generally, the two TADR data sets correlate well with the final lava morphology (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 correlation is even stronger. In the later stages, rubbly facies was emplaced under neither high TADR nor high variance conditions; however, qualitative 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 products, with both highly disrupted and continuous surfaces that exhibit inflation features (e.g., Keszthelyi et al., 2004; Vaucher et al., 2009; Jaeger et al., 2010; Voigt and Hamilton, 2018). The spatial 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 emplacement of spiny lava as local effusion rates decreased to become steadier and more continuous. Thus, critical information about effusion rate can be revealed for past fissure-fed eruptions by carefully examining the lava morphologies preserved within the geologic record.

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 Fulbright–Iceland. We thank Diego Coppola for providing the time-averaged discharge rate (TADR) data, and Ulrich Münzer and Iceland subglacial volcanoes interdisciplinary early warning system (IsViews) for the TerraSAR-X, TanDEM-X data, and UltraCam data.