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

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 km3 covering an area of ~84 km2. 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 (Th) 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 (Φrad) 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 m2) from a lava breakout has a close correspondence with a thermal camera measurement (1047 ◦C; occupying area of 4.52 m2). 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.


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 2 , with a bulk volume of 1.44 km 3  [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.
Remote Sens. 2017, 9, 151 2 of 24 [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.
Figure 1.The new Holuhraun lava field map by Icelandic Meteorological Office (after modification) [4].The Holuhraun lava field is situated on the flood plain, which is situated south of Askja volcano and north of Dyngjujökull, which is an outlet glacier from the Vatnajökull ice cap [1].

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 3 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 3 .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 3 /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 3 /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] Figure 1.The new Holuhraun lava field map by Icelandic Meteorological Office (after modification) [4].The Holuhraun lava field is situated on the flood plain, which is situated south of Askja volcano and north of Dyngjujökull, which is an outlet glacier from the Vatnajökull ice cap [1].

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 3 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 3 .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 3 /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 3 /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 3 /s [1,6].During this time, a lava pond <1 km 2 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 3 /s.In this phase, the lava transport was confined to closed lava pathways within flows.Over 19 km 2 of the flow field was resurfaced via surface breakouts from the closed pathways [1].
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 3 /s [1,6].During this time, a lava pond <1 km 2 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 3 /s.In this phase, the lava transport was confined to closed lava pathways within flows.Over 19 km 2 of the flow field was resurfaced via surface breakouts from the closed pathways [1].

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

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: where R x and R y are the radiances in bands x and y, respectively, (Wm −2 sr −1 m −1 ) after adjusting for atmospheric effects and surface emissivity; p is the pixel portion occupied by the hot component; R(λ x , T h ) and R λ y , T c are the radiances (Wm −2 sr −1 m −1 ) emitted for wavelengths λ x and λ 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]. 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.

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), 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.

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: where R SWIR and R TIR are the corrected spectral radiances at wavelength λ SWIR and λ TIR , respectively, L(λ SWIR ) and L(λ TIR ) are the spectral radiances at the sensor, L U (λ TIR ) is the atmospheric upwelling radiance, L R (λ SWIR ) is the atmospheric reflected radiance, τ is atmospheric transmissivity, and ε is the surface emissivity.In this case, we set the emissivity at 0.97 for the Holuhraun basaltic lava.R SWIR and R TIR will be used in Sections 4.1 and 4.2 both as an input for the calculation of TEI and for the dual-band method.

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 (R SWIR ) on the crust will be less than in the TIR spectral radiance (R TIR ) and vice versa on the active lava (R SWIR > R 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 where R SWIR and R TIR are the pixel corrected spectral radiances detected in the band 6 and band 10, respectively and R SWIR 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 Sections 5.1 and 5.3.
where R and R are the pixel corrected spectral radiances detected in the band 6 and band 10, respectively and R 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 Sections 5.1 and 5.

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, Th, Tc and p.In term of this study, we consider that Th relates to the temperature of recently active lava and recently crusted lava, while Tc 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].

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].where R and R are the pixel corrected spectral radiances detected in the band 6 and band 10, respectively and R 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 Sections 5.1 and 5.

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, Th, Tc and p.In term of this study, we consider that Th relates to the temperature of recently active lava and recently crusted lava, while Tc 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

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 (Φ rad ) for each pixel that contains lava can be estimated as where Φ rad is the radiant flux (W), σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W m −2 K −4 ), and A is the Landsat 8 pixel area, which is 900 m 2 .In this approach, we use the effective temperature model (T e ), expressed as 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, 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.

Convective Flux Estimation
Convective flux (Φ conv ) is estimated for the entire lava-flow field with a similar approach as above for Φ rad , i.e., given surface temperature T e and surfaces roughness H, Φ conv is calculated using the free convection formula given by [11,14] which is given by where the unit of Φ conv is W, h c = 5 W m −2 K −1 is the heat transfer coefficient for free convection [11], and T a is the ambient air temperature that is unaffected by eruption processes.In this work we use T a = 25 • C.

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 Then, re-arranged for ∆h: where ∆h is the crust thickness (m), M rad and M conv are radiative and convective flux densities (W m −2 ), respectively, k is the thermal conductivity, where we use 2.5 Wm −1 K −1 [14,31], and T 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].M rad and M conv are obtained by dividing Φ rad and Φ 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].

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 Then, re-arranged for ∆h: where ∆h is the crust thickness (m), M and M are radiative and convective flux densities (W m −2 ), respectively, k is the thermal conductivity, where we use 2.5 Wm −1 K −1 [14,31], and T 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].M and M are obtained by dividing Φ and Φ 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].

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

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.
Remote Sens. 2017, 9, 151 9 of 24 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

Spatial Distribution of Th and p
Table 3 shows the hot component temperature (T ) 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 Th 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 T , and p obtained from the Landsat 8 time series.The highest Th and lowest p mostly relate to active lava locations such as channels, lava ponds and breakouts.Meanwhile, the lowest Th 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.53m 2 , can be occupied by small freshly exposed lava patches or breakouts; this will be discussed further in Section 6.

Spatial Distribution of Th and p
Table 3 shows the hot component temperature (T 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 T 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.53m 2 , can be occupied by small freshly exposed lava patches or breakouts; this will be discussed further in Section 6.

Trend Th vs. p and TEI
We plotted Th as a function of p and TEI, respectively for 6 September 2014, as shown in Figure 8a,b, respectively.The scatter plot of Th vs. p shows a logarithmic decrease in Th 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 Th 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.51and 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 Th between 901 and 1185 °C and p between 0.10% and 1.30%.For sub-trend y, Th ranges between 400 and 900 °C while p in the range 0.10-7.5% and for sub-trend z, Th ranges between 346 and 750 °C while p remains in the range 0.10-13%.These trends are also apparent in scatter plots of Th versus TEI.The trends show large variation in Th within the subtrends.Examination of the scatter distribution for sub-trend y exhibits wide range of Th (346-1150 °C) being associated with a range 0.10-0.21 of TEI.For sub-trend z, Th ranges between 400 and 1160 °C while TEI is in the range 0.21-0.51.Meanwhile for sub-trend x, Th ranges between 901 and 1185 °C while exhibiting

Trend T h 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.51and 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.
Remote Sens. 2017, 9, 151 13 of 24 narrow range of TEI between 0.51 and 0.53.One can note from these results that maximum Th 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.

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

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).
Remote Sens. 2017, 9, 151 14 of 24 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).

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 groundbased 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.

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.

Discussions
In Sections 5.1 to 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 Tc 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

Discussions
In Section 5.1 , Sections 5.2 and 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 2 which is 0.33% of the pixel size, meanwhile FLIR shows a maximum temperature around 1047 • C for an area of 4.52 m 2 .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.
Remote Sens. 2017, 9, 151 16 of 24 measurement yields Th of 1096 °C in an area of 3.05 m 2 which is 0.33% of the pixel size, meanwhile FLIR shows a maximum temperature around 1047 °C for an area of 4.52 m 2 .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 Tc on Th for different thermal domains in this area.Figure 12 shows a logarithmic increase in Th as Tc increases [28]; therefore, Tc is an important parameter for deriving precise temperature estimates from the dual-band method.
(a) (b)  In Sections 5.4 and 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), measurement yields Th of 1096 °C in an area of 3.05 m 2 which is 0.33% of the pixel size, meanwhile FLIR shows a maximum temperature around 1047 °C for an area of 4.52 m 2 .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 Tc on Th for different thermal domains in this area.Figure 12 shows a logarithmic increase in Th as Tc increases [28]; therefore, Tc is an important parameter for deriving precise temperature estimates from the dual-band method.
(a) (b)  In Sections 5.4 and 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), In Sections 5.4 and 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.
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.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

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 2 ) from a lava breakout has a close correspondence with a thermal camera measurement (1047 °C; occupying area of 4.52 m 2 ).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.

Figure 3 .
Figure 3. R and R from Landsat 8 of the lava flow at Holuhraun, Iceland, on 6 September 2014.(a) Spectral radiances plot of SWIR (Band 6 and 7)-TIR (Band 10 and 11) and thermal mix model (hybrid Planck curve for the mixed pixel) (b) the band 6 detects the emissions from active lava (yellow-orange color); band 10 detects the emissions both from active lava and crust (red-green color).

Figure 4 .
Figure 4.The pixel mixture of the cool and the hot components.The region defined by the hot component temperature (Th) has pixel portion (p) and the region defined by the cool component temperature (Tc) has pixel portion (1-p): (a) A picture taken by Armann Hoskuldsson in Holuhraun, in 13 September 2014 showing the hot and the cool components; (b) The model of Landsat 8 mixed thermal pixel with 30 m pixel resolution.

Figure 3 .
Figure 3. R SWIR and R TIR from Landsat 8 of the lava flow at Holuhraun, Iceland, on 6 September 2014.(a) Spectral radiances plot of SWIR (Band 6 and 7)-TIR (Band 10 and 11) and thermal mix model (hybrid Planck curve for the mixed pixel) (b) the band 6 detects the emissions from active lava (yellow-orange color); band 10 detects the emissions both from active lava and crust (red-green color).

Figure 3 .
Figure 3. R and R from Landsat 8 of the lava flow at Holuhraun, Iceland, on 6 September 2014.(a) Spectral radiances plot of SWIR (Band 6 and 7)-TIR (Band 10 and 11) and thermal mix model (hybrid Planck curve for the mixed pixel) (b) the band 6 detects the emissions from active lava (yellow-orange color); band 10 detects the emissions both from active lava and crust (red-green color).

Figure 4 .
Figure 4.The pixel mixture of the cool and the hot components.The region defined by the hot component temperature (Th) has pixel portion (p) and the region defined by the cool component temperature (Tc) has pixel portion (1-p): (a) A picture taken by Armann Hoskuldsson in Holuhraun, in 13 September 2014 showing the hot and the cool components; (b) The model of Landsat 8 mixed thermal pixel with 30 m pixel resolution.

Figure 4 .
Figure 4.The pixel mixture of the cool and the hot components.The region defined by the hot component temperature (T h ) has pixel portion (p) and the region defined by the cool component temperature (T c ) has pixel portion (1-p): (a) A picture taken by Armann Hoskuldsson in Holuhraun, in 13 September 2014 showing the hot and the cool components; (b) The model of Landsat 8 mixed thermal pixel with 30 m pixel resolution.

Figure 5 .
Figure 5. Conductive model of the surface of lava and the basal thermal boundary layers for (a) an active smooth surface lava flow and (b) a crusted rough surface lava flow.(Adapted and modified from Harris [20] p. 222).

Figure 5 .
Figure 5. Conductive model of the surface of lava and the basal thermal boundary layers for (a) an active smooth surface lava flow and (b) a crusted rough surface lava flow.(Adapted and modified from Harris [20] p. 222).

Figure 6 .
Figure 6.(a) Spatial distribution of thermal eruption index (TEI) during the 6 September 2014 eruption: the four white dots show the location of the vents during the first phase of eruption [1]: (b) TEI differentiates four different thermal domains within the Holuhraun lava flow on 6 September 2014.TEI allows the active lava domain to be distinguished from crust domain and the non-volcanic hotspot domain.

Figure 6 .
Figure 6.(a) Spatial distribution of thermal eruption index (TEI) during the 6 September 2014 eruption: the four white dots show the location of the vents during the first phase of eruption [1]: (b) TEI differentiates four different thermal domains within the Holuhraun lava flow on 6 September 2014.TEI allows the active lava domain to be distinguished from crust domain and the non-volcanic hotspot domain.

Figure 7 .
Figure 7. Spatial distribution map of

Figure 8 .
Figure 8.(a) Scatter plot of Th vs. p; (b) Scatter plot of Th vs. TEI; (c) classification of lava thermal domain on 6 September 2014.

Figure 8 .
Figure 8.(a) Scatter plot of T h vs. p; (b) Scatter plot of T h vs. TEI; (c) classification of lava thermal domain on 6 September 2014.

Figure 9 .
Figure 9.Time series of the total radiant flux estimated from Landsat 8 for the 2014-2015 Holuhraun eruption.

Figure 9 .
Figure 9.Time series of the total radiant flux estimated from Landsat 8 for the 2014-2015 Holuhraun eruption.

Figure 10 .
Figure 10.(a) Crust thickness estimate derived from Landsat 8 by using the dual-band method during the eruption on 6 September 2014 (The black box shows the location of the field measurement; further details are given in Appendix D); (b) Thickness comparison between the satellite derived estimates and field measurements.

Figure 10 .
Figure 10.(a) Crust thickness estimate derived from Landsat 8 by using the dual-band method during the eruption on 6 September 2014 (The black box shows the location of the field measurement; further details are given in Appendix D); (b) Thickness comparison between the satellite derived estimates and field measurements.

Figure 11 .
Figure 11.(a) Lava breakout (yellow pixel) temperature from satellite from 2 December 2014, white star and arrow show the measurement location and the camera direction, respectively; (b) Lava breakout temperature measurement from a FLIR camera on 2 December 2014.The box shows the area of interest.

Figure 12 .
Figure 12.Th as function of Tc on 2 December 2014.Dots and solid lines are the solutions from dualband method and the logarithmic fitting respectively.

Figure 11 .
Figure 11.(a) Lava breakout (yellow pixel) temperature from satellite from 2 December 2014, white star and arrow show the measurement location and the camera direction, respectively; (b) Lava breakout temperature measurement from a FLIR camera on 2 December 2014.The box shows the area interest.

Figure 11 .
Figure 11.(a) Lava breakout (yellow pixel) temperature from satellite from 2 December 2014, white star and arrow show the measurement location and the camera direction, respectively; (b) Lava breakout temperature measurement from a FLIR camera on 2 December 2014.The box shows the area of interest.

Figure 12 .
Figure 12.Th as function of Tc on 2 December 2014.Dots and solid lines are the solutions from dualband method and the logarithmic fitting respectively.

Figure 12 .
Figure 12.T h as function of T c on 2 December 2014.Dots and solid lines are the solutions from dual-band method and the logarithmic fitting respectively.

Figure 13 .
Figure 13.Different values of H comparison for (a) radiant flux; (b) crust thickness, with the groundbased thickness in purple.

Figure 13 .
Figure 13.Different values of H comparison for (a) radiant flux; (b) crust thickness, with the ground-based thickness in purple.

Figure 13 .
Figure 13.Different values of H comparison for (a) radiant flux; (b) crust thickness, with the groundbased thickness in purple.

Figure A4 .
Figure A4.The Location of thickness measurement points, the arrows represent the order of points from 1 to 15.

Figure A4 .
Figure A4.The Location of thickness measurement points, the arrows represent the order of points from 1 to 15.

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

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

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

Table 3 .
Maximum, minimum and mean Th, p, and TEI values derived from the Landsat 8 time series.
DateTh min (°C) Th max (°C) Th average (°C) p Min p Max p Average TEI Min TEI Max

Table A1 .
Detail of the lava height measurement during 3-4 September 2014 and a comparison with satellite derived crust thickness measurement on 6 September 2014.

Table A1 .
Detail of the lava height measurement during 3-4 September 2014 and a comparison with satellite derived crust thickness measurement on 6 September 2014.