New Insights for Detecting and Deriving Thermal Properties of Lava Flow Using Infrared Satellite during 2014–2015 Effusive Eruption at Holuhraun, Iceland

Abstract

A new lava field was formed at Holuhraun in the Icelandic Highlands, north of Vatnajökull glacier, in 2014–2015. It was the largest effusive eruption in Iceland for 230 years, with an estimated lava bulk volume of ~1.44 km³ covering an area of ~84 km². Satellite-based remote sensing is commonly used as preliminary assessment of large scale eruptions since it is relatively efficient for collecting and processing the data. Landsat-8 infrared datasets were used in this study, and we used dual-band technique to determine the subpixel temperature (T_h) of the lava. We developed a new spectral index called the thermal eruption index (TEI) based on the shortwave infrared (SWIR) and thermal infrared (TIR) bands allowing us to differentiate thermal domain within the lava flow field. Lava surface roughness effects are accounted by using the Hurst coefficient (H) for deriving the radiant flux $\left({\mathsf{\Phi }}_{\mathrm{rad}}\right)$ and the crust thickness (Δh). Here, we compare the results derived from satellite images with field measurements. The result from 2 December 2014 shows that a temperature estimate (1096 °C; occupying area of 3.05 m²) from a lava breakout has a close correspondence with a thermal camera measurement (1047 °C; occupying area of 4.52 m²). We also found that the crust thickness estimate in the lava channel during 6 September 2014 (~3.4–7.7 m) compares closely with the lava height measurement from the field (~2.6–6.6 m); meanwhile, the total radiant flux peak is underestimated (~8 GW) compared to other studies (~25 GW), although the trend shows good agreement with both field observation and other studies. This study provides new insights for monitoring future effusive eruption using infrared satellite images.

1. Introduction

Holuhraun is a lava field in the Icelandic Highlands, north of Vatnajökull (Figure 1). The 2014–2015 Holuhraun lava field was created by basaltic fissure eruptions [1,2]. The eruptions lasted from 31 August 2014 to 27 February 2015 and formed a lava flow field covering 84 km², with a bulk volume of 1.44 km³ [1,2]. This eruption is the largest effusive eruption observed in the last 230 years in Iceland. Field observation during the 2014–2015 eruption has been documented in detail, i.e., Pedersen et al. [1] compiled a detailed evolution of the lava field and created a corresponding database containing information of the lava flows. Most of these field observations are related to the modes of lava transport and emplacement and thermal camera measurement, along with the mapping of the flow field growth and evolution [1,3]. This eruption offers an opportunity to improve our understanding of large effusive eruptions using satellite-based remote sensing. Here we present a new approach based on infrared satellite images to derive thermal properties within the lava field during eruption and then compare the results with field measurement.

Study Area

Iceland is one of the most active volcanic regions on Earth; volcanic eruptions occur, on average, every four to five years, and produce more than 5 km³ magma per century [5]. The 6-month long eruption at Holuhraun 2014–2015 was the largest effusive eruption in Iceland in 230 years, with an estimated bulk lava volume of about 1.44 km³. The 2014–2015 Holuhraun lava flows are emplaced on the floodplain 0–6 km from the Dyngjajökull glacier (Figure 1) occupying a relatively flat area. The eruption had an average discharge of about 77 m³/s, making it the longest effusive eruption observed in modern times with such a flux [1,6]. According to the study of Pedersen et al. [1], the eruption was divided into three phases based on the lava field evolution, as follows:

• Phase 1: Open channel lava pathways: 31 August to mid-October (Figure 2a);
• Phase 2: Lava pond formation: Mid-October to end-November (Figure 2b);
• Phase 3: Tube-fed lava pathways: early December to 27 February (Figure 2c).

The first phase of the 2014–2015 Holuhraun eruption was dominated by open lava channels. This phase had a discharge ~350–100 m³/s. The eruption began on a 1.8 km long fissure feeding up to 500 m wide, incandescent sheets of slabby pahoehoe [1,2]. Apart from fire fountaining in the first few weeks, in 15 days, the volcanic activity had formed an open channels lava flow that advanced 17–18 km towards NNE and subsequently the lava morphology changed to rubbly and aa types [1,7]. The second phase had a discharge ranging from 50 to 100 m³/s [1,6]. During this time, a lava pond <1 km² was formed at a distance of 0.8 km east of the vents [1]. This pond became the main point of lava distribution, controlling the emplacement of the lava flows [1,2]. Towards the end of this phase, the open channel in the first phase was inflated due to new lava injections into the previously active lava channel lifting the channel [1]. The final (third) phase from December to the end of February had a mean discharge <50 m³/s. In this phase, the lava transport was confined to closed lava pathways within flows. Over 19 km² of the flow field was resurfaced via surface breakouts from the closed pathways [1].

2. Infrared Remote Sensing in Volcanoes

Satellite-based infrared remote sensing data are increasingly being used to monitor active volcanoes around the world [8,9,10]. Monitoring volcanoes by infrared remote sensing is essential to improve understanding of active volcanoes. For example, thermal monitoring of volcanoes is necessary in order to understand the volcanic eruption processes. Effusive eruption activity can change rapidly over time as new lava flow form and develop. In many cases, thermal observations of active eruptions from both ground and aircraft are risky and difficult, especially when the lava covers a large area. Satellite-based remote sensing provides high temporal resolution infrared data that are suitable for monitoring long effusive eruption such as in the new lava field at Holuhraun, Iceland. Such a tool provides data in which effusive events are detectable and changes in the eruption style and evolution of activity can be identified regardless of the coarse spatial resolution of the data [11]. Satellite-based remote sensing can be used to monitor thermal activity within lava fields [9]. In the last 25 years, there have been several methods that have been used for deriving eruption activity information from infrared remotely-sensed data to estimate the thermal structures of hot volcanic features such as active lava flows [12,13,14,15,16], heat flux [11,17,18], effusive rate [10,11,17,18] and crust thickness [10,11]. In 1981, Dozier [19] developed a method involving a solution of simultaneous equations that allows the calculation of the ‘sub-pixel’ coverage and temperature of cool and hot components. This method is called the dual-band method, and involves two distinct infrared bands to formulate a system of two equations from the simultaneous solution of the Planck equation in each band as shown below:

$${\mathrm{R}}_{\mathrm{x}}=\mathrm{pR}\left({\mathsf{\lambda }}_{\mathrm{x}},{\mathrm{T}}_{\mathrm{h}}\right)+\left(1-\mathrm{p}\right)\mathrm{R}\left({\mathsf{\lambda }}_{\mathrm{x}},{\mathrm{T}}_{\mathrm{c}}\right)$$

(1)

$${\mathrm{R}}_{\mathrm{y}}=\mathrm{pR}\left({\mathsf{\lambda }}_{\mathrm{y}},{\mathrm{T}}_{\mathrm{h}}\right)+\left(1-\mathrm{p}\right)\mathrm{R}\left({\mathsf{\lambda }}_{\mathrm{y}},{\mathrm{T}}_{\mathrm{c}}\right)$$

(2)

where R_x and R_y are the radiances in bands x and y, respectively, (Wm⁻² sr⁻¹ m⁻¹) after adjusting for atmospheric effects and surface emissivity; p is the pixel portion occupied by the hot component; $\mathrm{R}\left({\mathsf{\lambda }}_{\mathrm{x}},{\mathrm{T}}_{\mathrm{h}}\right)$ and $\mathrm{R}\left({\mathsf{\lambda }}_{\mathrm{y}},{\mathrm{T}}_{\mathrm{c}}\right)$ are the radiances (Wm⁻² sr⁻¹ m⁻¹) emitted for wavelengths ${\mathsf{\lambda }}_{\mathrm{x}}$ and ${\mathsf{\lambda }}_{\mathrm{y}}$, at surface temperatures T_h (hot component) and T_c (cool component), respectively. The dual-band method can be applied if the two bands of the short-wave infrared (SWIR) and the thermal infrared (TIR) data are available [20], whereby any two of the unknowns, T_c, T_h and p, can be estimated if the third is assumed. This method has been successfully applied by several researchers [5,6,7,8,9,10,13,14]. Harris et al. and Lombardo et al. [11,16] used dual band method to retrieve the crust and the hot cracks temperature for active lava flows in Mt. Etna. They used band 5 (1.55–1.75 µm) and 7 (2.08–2.35 µm) from Landsat Thematic Mapper (TM) and assumed T_h to estimate T_c and p. In this study, we develop a new spectral index for Landsat 8, named the thermal eruption index (TEI), based on the SWIR and TIR bands (bands 6 and 10). The purpose of the TEI consists mainly of two parts: (1) as a threshold for differentiating between different thermal domains; and (2) applying a dual-band method to estimate subpixel temperature within thermal domains and differentiating between the types of lava surface. The active lava surface has thermal domain complexity and could contain more than one thermal component [13,21]; here we use two thermal component scenarios, with T_h as temperature of lava surface and T_c as temperature surrounding the lava for different thermal domains. The remainder of this manuscript is organized as follows. Section 3 describes the datasets and the preprocessing. The proposed methodology for TEI, the dual-band method, the estimation of radiant flux (Φ_rad), and the crust thickness model of lava flow (Δh) is explained in Section 4. We also discuss the effect of lava surface roughness using Hurst coefficient (H) on Φ_rad and Δh. Further results are arranged in Section 5, and in Section 6 we give a brief discussion. In Section 7, we give a summary of our work.

3. Datasets and Preprocessing

3.1. Datasets

Remote sensing observations were made using Landsat 8 Level 1 product band 6 (1.56–1.66 µm) and band 10 (10.60–11.19 µm). The eruption was well monitored from Landsat 8, and although Landsat 8 only has a temporal resolution of once every 16 days (making it of limited value for making time series studies of the eruption), the spatial resolution makes it a good tool to derive thermal properties within the lava flow. Landsat 8 has different spatial resolution for band 6 and band 10: band 6 has 30 m spatial resolution and band 10 has 100 m spatial resolution (resampled to 30 m). The selection of band 6 and 10 are considered to minimize oversaturation effects over active lava flows. According to Blackett [22], Landsat 8 has an enhanced dynamic range compared to Landsat ETM+: this means that temperatures of up to 747.9 K can be detected without saturation in band 6 as compared with those for the corresponding ETM+ band 724.5 K. Acquisition dates are selected according to the availability and quality of data covering the eruption (

Table 1. Product ID and the dates of the Landsat 8 datasets that were used in this study.

Product ID	Date
LC82170152014249LGN00	6 September 2014
LC82180142014272LGN00	29 September 2014
LC82180142014288LGN00	15 October 2014
LC80642292014297LGN00	24 October 2014
LC80642302014329LGN00	25 November 2014
LC80652292014336LGN00	2 December 2014
LC80652292014352LGN00	18 December 2014
LC80652292015003LGN00	3 January 2015
LC82180142015035LGN00	4 February 2015

), we only took the data where cloud coverage is minimal. The data can be downloaded from the U.S Geological Survey (USGS) website (https://earthexplorer.usgs.gov). In our work, we subset Landsat 8 images into 562 by 333 and then converted the satellite-recorded digital numbers (DN) to sensor radiance for both SWIR and TIR bands. In this study, we use radar Sentinel 1A data from 18 October 2014 to derive Hurst coefficient from the lava. This data represents roughness of the lava and can be downloaded from the website (https://scihub.copernicus.eu/dhus/#/home). We also use thermal camera (FLIR) measurement during 2 December 2014 that overlaps with the satellite data and theodolite lava height measurement in the field for result comparison.

3.2. Atmospheric and Emissivity Correction

The MODTRAN model atmosphere is used for atmospheric correction in this study [13,23]. Since emissivity and atmospheric effects will vary by wavelength [20], corrections are needed for both SWIR and TIR as follows:

$${\mathrm{R}}_{\mathrm{SWIR}}=\left[\mathrm{L}\left({\mathsf{\lambda }}_{\mathrm{SWIR}}\right)-{\mathrm{L}}_{\mathrm{R}}\left({\mathsf{\lambda }}_{\mathrm{SWIR}}\right)\right]/\mathsf{\tau }\mathsf{\epsilon }$$

(3)

$${\mathrm{R}}_{\mathrm{TIR}}=\left[\mathrm{L}\left({\mathsf{\lambda }}_{\mathrm{TIR}}\right)-{\mathrm{L}}_{\mathrm{U}}\left({\mathsf{\lambda }}_{\mathrm{TIR}}\right)\right]/\mathsf{\tau }\mathsf{\epsilon }$$

(4)

where ${\mathrm{R}}_{\mathrm{SWIR}}$ and ${\mathrm{R}}_{\mathrm{TIR}}$ are the corrected spectral radiances at wavelength ${\mathsf{\lambda }}_{\mathrm{SWIR}}$ and ${\mathsf{\lambda }}_{\mathrm{TIR}}$, respectively, $\mathrm{L}\left({\mathsf{\lambda }}_{\mathrm{SWIR}}\right)$ and $\mathrm{L}\left({\mathsf{\lambda }}_{\mathrm{TIR}}\right)$ are the spectral radiances at the sensor, ${\mathrm{L}}_{\mathrm{U}}\left({\mathsf{\lambda }}_{\mathrm{TIR}}\right)$ is the atmospheric upwelling radiance, ${\mathrm{L}}_{\mathrm{R}}\left({\mathsf{\lambda }}_{\mathrm{SWIR}}\right)$ is the atmospheric reflected radiance, $\mathsf{\tau }$ is atmospheric transmissivity, and $\mathsf{\epsilon }$ is the surface emissivity. In this case, we set the emissivity at 0.97 for the Holuhraun basaltic lava. ${\mathrm{R}}_{\mathrm{SWIR}}$ and ${\mathrm{R}}_{\mathrm{TIR}}$ will be used in Section 4.1 and Section 4.2 both as an input for the calculation of TEI and for the dual-band method.

4. Method

4.1. Thermal Eruption Index (TEI)

In this study, TEI is developed by using the SWIR and the TIR bands from the medium spatial resolution satellite Landsat 8. The method uses the sensitivity difference between SWIR (band 6) and TIR (band 10) to detect pixel hot spots instead of using mid-infrared (MIR) as in NTI (Normalized Thermal Index) [22,24]. The objective of TEI is to provide a new variant for a hotspot thermal index derived by using data from the medium spatial resolution satellite Landsat 8. We derive TEI based on image observation and define the empirical formula for this observation. TEI is based on the principle that the SWIR spectral radiance (${\mathrm{R}}_{\mathrm{SWIR}}$) on the crust will be less than in the TIR spectral radiance (${\mathrm{R}}_{\mathrm{TIR}}$) and vice versa on the active lava (${\mathrm{R}}_{\mathrm{SWIR}}$ > ${\mathrm{R}}_{\mathrm{TIR}}$) as shown in Figure 3a. Figure 3b shows band 6 and band 10 of Landsat 8 from the 6 September 2014 Holuhraun eruption, the active lava pixels are emitting more spectral radiance in both band 6 and band 10; meanwhile the crust pixels are emitting more spectral radiance only in TIR. Therefore TEI has higher values in the active lava than in crust. This index uses the square of the TIR spectral radiance and the maximum of the SWIR spectral radiance to differentiate between the thermal domains. TEI is expressed as

$$\mathrm{TEI}=\frac{{\mathrm{R}}_{\mathrm{SWIR}}-\frac{{\left({\mathrm{R}}_{\mathrm{TIR}}\right)}^2}{10{\mathrm{R}}_{\mathrm{SWIR}\mathrm{MAX}}}}{{\mathrm{R}}_{\mathrm{SWIR}}+\frac{{\left({\mathrm{R}}_{\mathrm{TIR}}\right)}^2}{10{\mathrm{R}}_{\mathrm{SWIR}\mathrm{MAX}}}}\frac{{\left({\mathrm{R}}_{\mathrm{TIR}}\right)}^2}{{(\frac{{\mathrm{R}}_{\mathrm{SWIR}\mathrm{MAX}}}{3})}^2}$$

(5)

where ${\mathrm{R}}_{\mathrm{SWIR}}$ and ${\mathrm{R}}_{\mathrm{TIR}}$ are the pixel corrected spectral radiances detected in the band 6 and band 10, respectively and ${\mathrm{R}}_{\mathrm{SWIR}\mathrm{MAX}}$ are the maximum spectral radiances detected in band 6 for each scene. In this study, we applied the dual band method to automatically calculate the hot component temperature within the region defined by the hotspot threshold (TEI > 0.10). This selected value is explained in more detail in Section 5.1 and Section 5.3.

4.2. Dual-Band Method

The dual-band method relies on a system of two equations (Equations (1) and (2)) and requires the assumption of one of three unknowns, T_h, T_c and p. In term of this study, we consider that T_h relates to the temperature of recently active lava and recently crusted lava, while T_c relates to the surrounding temperature that is influenced by the eruption processes (Figure 4). both of these can be assumed; meanwhile, p is typically difficult to assume accurately and therefore has to be derived [22].

In this study, we set T_c equal to the lowest brightness temperature detected in TIR for each thermal domain considered, with T_c = 25 °C in the surrounding warm crust, T_c = 50 °C in the surrounding hot crust, and T_c = 85 °C in the active lava. These assumptions will be suitable in situations where different thermal domains (active lava and crust) within the lava flow are clearly separable. We solve p by iterating on T_h, until p(R_SWIR) = p(R_TIR) then we can rearrange Equations (1) and (2) into

$$\mathrm{p}\left({\mathrm{R}}_{\mathrm{SWIR}}\right)=\frac{{\mathrm{R}}_{\mathrm{SWIR}}-\mathrm{R}\left({\mathsf{\lambda }}_{\mathrm{SWIR}},{\mathrm{T}}_{\mathrm{c}}\right)}{\mathrm{R}\left({\mathsf{\lambda }}_{\mathrm{SWIR}},{\mathrm{T}}_{\mathrm{h}}\right)-\mathrm{R}\left({\mathsf{\lambda }}_{\mathrm{SWIR}},{\mathrm{T}}_{\mathrm{c}}\right)}$$

(6)

$$\mathrm{p}\left({\mathrm{R}}_{\mathrm{TIR}}\right)=\frac{{\mathrm{R}}_{\mathrm{TIR}}-\mathrm{R}\left({\mathsf{\lambda }}_{\mathrm{TIR}},{\mathrm{T}}_{\mathrm{c}}\right)}{\mathrm{R}\left({\mathsf{\lambda }}_{\mathrm{TIR}},{\mathrm{T}}_{\mathrm{h}}\right)-\mathrm{R}\left({\mathsf{\lambda }}_{\mathrm{TIR}},{\mathrm{T}}_{\mathrm{c}}\right)}$$

(7)

4.3. Radiant Flux Estimation

Radiation is the most direct heat flux to estimate. For rough lava surface (aa and brecciated surface), not all the radiation can escape from the lava surface because of surface scattering. Therefore, in this paper, we propose to use the Hurst coefficient (H) [25,26] to describe the surface roughness of lava, so that the actual radiation emitted is reduced due to the fractal model. Following this model, the radiant flux $\left({\mathsf{\Phi }}_{\mathrm{rad}}\right)$ for each pixel that contains lava can be estimated as

$${\mathsf{\Phi }}_{\mathrm{rad}}={\mathsf{\epsilon }\mathsf{\sigma }\mathrm{HAT}}_{\mathrm{e}}{}^4$$

(8)

where ${\mathsf{\Phi }}_{\mathrm{rad}}$ is the radiant flux (W), $\mathsf{\sigma }$ is the Stefan–Boltzmann constant (5.67 × 10⁻⁸ W m⁻² K⁻⁴), and A is the Landsat 8 pixel area, which is 900 m². In this approach, we use the effective temperature model (T_e), expressed as

$${\mathrm{T}}_{\mathrm{e}}={\left({\mathrm{pT}}_{\mathrm{h}}{}^4+\left(1-\mathrm{p}\right){\mathrm{T}}_{\mathrm{c}}{}^4\right)}^{1/4}$$

(9)

which is the average surface temperature of lava for the two thermal components present on the lava flow surface [27,28]. H is the Hurst coefficient (0 < H < 1), where a higher H means a smoother surface. We set H for each thermal domain within lava field as shown in

Table 2. The values of the Hurst coefficient for different thermal domains.

H’	H (after Normalized)	Description	Thermal Domain
0.43	0.21	Very rough surface (Aa flow, large, tilted, spinny pates [30]).	Warm crust
0.70	0.35	Rough surface (Aa flow, small spinny plates [30]).	Hot crust
0.89	0.44	Smooth surface (Sheet pahoehoe flow, channel, roppy structure [30]).	Active lava

, these values were derived from radar image speckle pattern transect (H’) [29] and we normalized with 0.5 according the study by Shepard et al. [30]. Further details about the derivation of H from radar image are explained in Appendix A.

4.4. Convective Flux Estimation

Convective flux $({\mathsf{\Phi }}_{\mathrm{conv}})$ is estimated for the entire lava-flow field with a similar approach as above for ${\mathsf{\Phi }}_{\mathrm{rad}}$, i.e., given surface temperature T_e and surfaces roughness H, ${\mathsf{\Phi }}_{\mathrm{conv}}$ is calculated using the free convection formula given by [11,14] which is given by

$${\mathsf{\Phi }}_{\mathrm{conv}}={\mathrm{Ah}}_{\mathrm{c}}\mathrm{H}\left({\mathrm{T}}_{\mathrm{e}}-{\mathrm{T}}_{\mathrm{a}}\right)$$

(10)

where the unit of ${\mathsf{\Phi }}_{\mathrm{conv}}$ is W, ${\mathrm{h}}_{\mathrm{c}}=5$ W m⁻² K⁻¹ is the heat transfer coefficient for free convection [11], and ${\mathrm{T}}_{\mathrm{a}}$ is the ambient air temperature that is unaffected by eruption processes. In this work we use ${\mathrm{T}}_{\mathrm{a}}$ = 25 °C.

4.5. Crust Thickness Model

The crust thickness Δh is calculated by assuming that the conductive flux density across the surface crust is equal to the total of the radiative and convective flux densities leaving the same surface of lava [20,23], so that

$${\mathrm{M}}_{\mathrm{rad}}+{\mathrm{M}}_{\mathrm{conv}}=-\mathrm{k}\frac{{\mathrm{T}}_{\mathrm{i}}-{\mathrm{T}}_{\mathrm{e}}}{\Delta \mathrm{h}}$$

(11)

Then, re-arranged for $\Delta \mathrm{h}$:

$$\Delta \mathrm{h}=-\mathrm{k}\frac{\Delta \mathrm{T}}{{\mathrm{M}}_{\mathrm{rad}}+{\mathrm{M}}_{\mathrm{conv}}}$$

(12)

where $\Delta \mathrm{h}$ is the crust thickness (m), ${\mathrm{M}}_{\mathrm{rad}}$ and ${\mathrm{M}}_{\mathrm{conv}}$ are radiative and convective flux densities (W m⁻²), respectively, $\mathrm{k}$ is the thermal conductivity, where we use 2.5 Wm⁻¹ K⁻¹ [14,31], and ${\mathrm{T}}_{\mathrm{i}}$ is the temperature of the lava flow interior. In this study, we use an interior temperature of 1128 °C (for lava outside vent) and 1200 °C (for lava surrounding vent); these values were selected according to thermocouple measurements for freshly exposed patches of lava in Holuhraun on the 19 and 20 November 2014 [3]. ${\mathrm{M}}_{\mathrm{rad}}$ and ${\mathrm{M}}_{\mathrm{conv}}$ are obtained by dividing ${\mathsf{\Phi }}_{\mathrm{rad}}$ and ${\mathsf{\Phi }}_{\mathrm{conv}}$ by the pixel area A. Figure 5 depicts the physical meaning of Δh [10,20]. On rough surface (aa and brecciated surface) Δh is the thickness of the thermal boundary beneath a thermally mixed surface rubble [10,20]. Meanwhile, since there is no surface rubble in smooth surface (brittle layer and thin), we can assume that Δh is the complete thickness of the thermal boundary of the lava surface [10,20].

5. Results

5.1. TEI Hotspots Anomaly

We detected different types of hotspots associated with the 2014–2015 Holuhraun eruption by using TEI. In this section, we review some of these results to illustrate the utility of the TEI for differentiating thermal anomalies within the lava field. As noted in Section 4.1, TEI detects hotspots with TEI > 0.10, but it does not discriminate between the different types of the hotspot domains that occurs in the lava field. As a result, the interpretation can be difficult, as the active lava, hot crust and warm crust could be mistaken for eruptive activity. Therefore, accurate lava evolution where the hotspots occur becomes critical for determining the domain of TEI anomaly. Thus, we selected the image from 6 September 2014 to assess the TEI value, since the image has a clear difference between active lava and recently emplaced crust. Figure 6a shows the spatial distribution of TEI on 6 September 2014. In general, TEI values detect hotspots in the range from 0.10 to 0.53. We distinguish the two main thermal domains within the lava flow field. The first is the active lava domain, which is characterized by high TEI. The second is the crust domain surrounding active lava, characterized by TEI below active lava domain but exceeding the hotspot threshold (>0.10). Figure 6b shows that the active lava domains have TEI value of >0.51: this value is related to high emitted radiance in SWIR and TIR, whilst crust domains have value ranging from 0.10 to 0.51. The crust zone itself is divided into two sub-domains: (1) hot crust domain, having values ranging from 0.21 to 0.51; (2) warm crust domain, having values ranging from 0.10 to 0.21, This also can be seen during 29 September and 24 October (See Appendix B). Clearly, TEI allows better discrimination within the lava flow; active lava, crust and non-volcanic hotspot, offering the possibility that different mode of lava flows can be automatically discriminated based of TEI threshold as will be discussed in Section 5.3.

5.2. Spatial Distribution of Th and p

Table 3. Maximum, minimum and mean Th, p, and TEI values derived from the Landsat 8 time series.

Date	Th min (°C)	Th max (°C)	Th average (°C)	p Min	p Max	p Average	TEI Min	TEI Max
6 September 2014	346	1185	685	0.14%	13%	2.2%	0.102	0.53
29 September 2014	442	1199	885	0.10%	6%	1%	0.100	0.53
15 October 2014	526	1195	816	0.11%	5%	1.2%	0.101	0.53
24 October 2014	347	1195	705	0.13%	13%	2%	0.100	0.53
25 November 2014	392	1201	901	0.11%	7%	0.08%	0.102	0.53
2 December 2014	368	1208	770	0.10%	10%	1.4%	0.101	0.53
18 December 2014	358	1191	677	0.17%	10%	2.1%	0.100	0.53
3 January 2015	344	1176	689	0.15%	13%	2%	0.102	0.53
4 February 2015	487	1154	790	0.12%	4%	1%	0.101	0.53

shows the hot component temperature (${\mathrm{T}}_{\mathrm{h}}$) and the pixel portion of the hot component (p) computed by using the TEI based dual band-method for nine Landsat 8 time series during the eruption. Our solutions of T_h ranges between 344 °C and 1208 °C, and p ranges between 0.10 and 13% with an average of 769 °C and 1.5%, respectively. Figure 7 shows the spatial distribution map of ${\mathrm{T}}_{\mathrm{h}}$, and p obtained from the Landsat 8 time series. The highest T_h and lowest p mostly relate to active lava locations such as channels, lava ponds and breakouts. Meanwhile, the lowest T_h and highest p corresponds to crust zone in the edge of active lava and flow fronts that indicates lava cooling and formed thick crust, as lava cools the crust component increases and progressively thickens from the vent towards the distal end. Such results are broadly consistent with observed emplacement processes of aa flows on the 6 September 2014. The result points out that lowest p (0.10–0.17%) of a pixel), equivalent to 0.90–1.53 m², can be occupied by small freshly exposed lava patches or breakouts; this will be discussed further in Section 6.

5.3. Trend Th vs. p and TEI

We plotted T_h as a function of p and TEI, respectively for 6 September 2014, as shown in Figure 8a,b, respectively. The scatter plot of T_h vs. p shows a logarithmic decrease in T_h as p increases: this means that the higher the temperature, the smaller the pixel portion. This general trend is in good agreement with theoretical models [13,14,16,21]. There are three sub-trends (x, y and z) identified on the scatter plots (Figure 8a,b), and these sub trends allow for classification of lava thermal domain as shown in Figure 8c. The same trends are observed on open channel flow during 29 September and 24 October (See Appendix C). This classification shows that the active lava domain is associated with T_h in the range 901–1208 °C and p in the range 0.1–0.13% with TEI > 0.51, while the hot crust domain is in the range of 400–900 °C and p in the range of 0.20–7.5% with TEI 0.21–0.51 and warm crust domain is in the range of 346–750 °C and p is in the range of 0.30–13% with TEI 0.10–0.22. This classification result has good agreement with lava flow field observations. In Figure 8a, we observe a nearly linear decrease for sub-trend x; on the other hand, the scatter plot associated with sub-trend y and z show logarithmic decrease as in the general trend. The scatter distribution for sub-trend x has its support in region enclosed by T_h between 901 and 1185 °C and p between 0.10% and 1.30%. For sub-trend y, T_h ranges between 400 and 900 °C while p in the range 0.10–7.5% and for sub-trend z, T_h ranges between 346 and 750 °C while p remains in the range 0.10–13%. These trends are also apparent in scatter plots of T_h versus TEI. The trends show large variation in T_h within the subtrends. Examination of the scatter distribution for sub-trend y exhibits wide range of T_h (346–1150 °C) being associated with a range 0.10–0.21 of TEI. For sub-trend z, T_h ranges between 400 and 1160 °C while TEI is in the range 0.21–0.51. Meanwhile for sub-trend x, T_h ranges between 901 and 1185 °C while exhibiting narrow range of TEI between 0.51 and 0.53. One can note from these results that maximum T_h could occur in TEI < 0.51. This can happen for several reasons: (a) The effect of plumes from the eruption that mix within the lava pixels. This means that, due to the sensibility of band 6 [32], the plume was detected as high radiance in band 6. This makes TEI > 0.1 (hotspot threshold) and the dual-band produces very high temperatures, while TEI < 0.51: (b) A different spatial resolution between band 6 (30 m) and band 10 (100 m) causes lower TEI value due to the small fragments of high emission lava in band 6 that is not detected in band 10.

5.4. Radiant Flux Time Series during 2014–2015 Holuhraun Eruption

Figure 9 depicts the temporal evolution of the radiant flux (total energy radiated) for nine observations of the 2014–2015 Holuhraun effusive eruptions. As mentioned in Section 3.1, Landsat 8 has limited value for making time series studies of the eruption, since the temporal resolution is low and there are quite a lot of flux gaps. However, the trend shows good agreement with results from MODIS done by Wright et al. [33] during the first 101 days of eruption, although our total radiant flux peak is underestimated compared to Wright et al. [33] since we consider the effect of the Hurst coefficient (H) parameter that affects the total flux value (discussed in Section 6). The maximum peak detected is on day 7 of the eruption (6 September 2014), with a total flux of 7.8 GW, and then it continues to decline until it reaches a low level of 1.6 GW on day 46 before it rises again to a new peak of 3.2 GW on day 55. It then decreases to a new low on day 87 with 1.3 GW and then it increases to a new peak on day 110 with 1.6 GW. Then, it continues to decrease until day 158 with 0.2 GW, indicating that the eruption has almost stopped. This total flux peak trend is also in good agreement with field observations [1,2]: the first phase of eruption is related to the evolution of the vent system from a fissure with discrete vents distributed along the length of the fissure (day 7) decreasing in number with time to eruption from a single source vent (day 30) and formed a lava pond in the second phase (day 46) (Figure 7a–c). Then, on day 55, the pond becomes the main lava distributor and the pond size continues to decrease [1,2] until the final phase of eruption starts on day 110 (Figure 7d–f). In the final phase of eruption, during days 110, 126 and 158, the flow field was mostly formed via surface breakouts [1,2] until the eruption stopped (Figure 7g–i).

5.5. Crust Thickness Model of Lava Flow

Figure 10a shows the crust thickness estimate of lava flow for 6 September 2014. The thickness estimates range from 2 to15 m. The thickest crust is located along the edge of lava flow and in lava flow front. Meanwhile the thinnest part is mostly located in the active lava flow (lava channel and breakouts). This result also has a good agreement with Rossi et al. [34] that shows using TanDEM-X data, that the thickness in the lava field ranges between 0 and15 m during 9 September 2014 with the thicker part along the edge and the thinner parts in the middle of the channel. For comparison, 15 points of lava height measurement were done using Theodolite during 3 to 4 September 2014 (location measurement points are in Appendix D). Figure 10b shows the comparison between thickness measurement from the field and the satellite derived estimate. We found that the ground-based thickness measurement are closer to the satellite derived estimates in the middle of lava channel. Meanwhile, along the edge of lava flow, the satellite derived thickness estimates are thicker than the ground-based measurements. Presumably, this is because the lava in the edge cools down and fully develops into crust on satellite, since there are temporal gaps between the field observations and the satellite derived estimates. On the other hand, the lava in the middle of the channel is still active and has not completely cooled, and this leads to small thickness difference between the field measurement and the satellite derived estimates. We also note that we only derive the crust thickness and not complete thickness of the lava, since there are still have fluid interior layer beneath the crust, especially for active lava flow, as seen in Figure 5.

6. Discussions

In Section 5.1 , Section 5.2 and Section 5.3, we have used TEI to automatically detect and derive thermal properties from infrared remotely-sensed data within the hotspot region defined by TEI > 0.10. This value provides encouragement that the TEI method yields robust estimates of hotspot anomalies during eruption. On the other hand, differentiating between thermal domains offers new possibilities to use different T_c setting for each domain to derive sub pixel temperature using the dual-band method. However, as explained by Lombardo and Buongiorno [13], the dual-band method provides a rough approximation to the thermal model when only two infrared bands are available. Figure 11a,b show comparison of active lava temperature obtained from satellite and forward looking infrared (FLIR) camera measurements during 2 December 2014 lava breakout in 64°54.644’N, 16°42.931’W. This comparison shows a good agreement for both satellite and field measurements. Satellite measurement yields T_h of 1096 °C in an area of 3.05 m² which is 0.33% of the pixel size, meanwhile FLIR shows a maximum temperature around 1047 °C for an area of 4.52 m². However, this comparison is a rough estimation, since the temperature is not uniform within the active lava. Therefore, we recommend further improvement by using more than two thermal components. We also performed an experiment to study the effect of different T_c on T_h for different thermal domains in this area. Figure 12 shows a logarithmic increase in T_h as T_c increases [28]; therefore, T_c is an important parameter for deriving precise temperature estimates from the dual-band method.

In Section 5.4 and Section 5.5, we found that both the radiant flux estimate and crust thickness estimate agree closely with the field observations trends. In Figure 13a,b, we conduct an experiment by varying H: for very rough lava (H = 0.3), rough lava (H = 0.7) and perfectly smooth lava (H = 1). Smoother lava will produce the highest total radiant flux, but on the other hand, the crust becomes thinner. The higher H (0.7–1) agrees with the radiant flux peak result of Wright et al. [33] (~25 GW), since the effect of surface roughness are not accounted in their study. Interestingly, in the channel and northern part of the crust (Figure A4, points 1 to 6) results are produced that are closer to the field measurements; this means that the lava in those thermal domains is mostly dominated by rough surface (lower H). This is in good agreement with field observations that show that brecciated surface is dominant in Holuhraun. Meanwhile, the southern part of the crust (Figure A4, points 10 to 15) is overestimated, which means that the H is higher (smoother surface) than we estimate. Simply, as we mentioned earlier, this is due to the temporal gap between field measurement and satellite. In Figure 14a,b, it is shown that the temperature of the channel decreases with distance from the vents and crust thickness increases from the vents on 6 September 2014 (Figure 7a); This trend has good agreement with past work from Oppenheimer during the Lonquimay eruption (Chile, 1989) [10] which showed such trends. However, the H parameter is still open for discussion, and we recommend performing an alternative method to determine H (e.g., airborne light detection and ranging (LIDAR) and terrestrial laser scanning) and varying H based on LIDAR model during future eruption in order to get a better estimate.

7. Conclusions

In this paper, we introduce a new spectral index called the thermal eruption index (TEI) based on the SWIR and TIR bands, allowing us to differentiate thermal domains within the lava flow field. TEI detects hotspots with TEI > 0.10: this value provides encouragement that the TEI method yields robust estimates of hotspot anomalies during eruption. Two main thermal domains were distinguished within the lava flow field. The first is the active lava domain, which is characterized by high TEI (0.51). The second is the crust domain surrounding active lava, characterized by TEI below active lava domain but exceeding the hotspot threshold (0.10–0.51). The result from 2 December 2014 shows that a temperature estimate (1096 °C; occupying area of 3.05 m²) from a lava breakout has a close correspondence with a thermal camera measurement (1047 °C; occupying area of 4.52 m²). This paper also considered effect of lava surface roughness effects by using the Hurst coefficient (H) for deriving the radiant flux (Φ_rad) and the crust thickness (Δh), where the higher H (smoother surface) produce thinner crust meanwhile the lower H (rough surface) will produce thicker crust. Crust thickness in the lava channel during 6 September 2014 (~3.4–7.7 m) compares closely with the lava height measurement from the field (~2.6–6.6 m); meanwhile, the total radiant flux peak is underestimated (~8 GW) compared to other studies (~25 GW), although the trend shows good agreement with both field observation and other studies. These results show that the proposed techniques were successfully applied to Landsat 8 on SWIR and TIR datasets from 2014–2015 Holuhraun eruptions. In future work, the proposed methods will be applied to other satellite/airborne datasets which have both SWIR and TIR band and consider alternative method to determined H (e.g., airborne LIDAR and terrestrial laser scanning). This study provides new insights for monitoring future effusive eruption using infrared satellite images.