Lava Flow Roughness on the 2014–2015 Lava Flow-Field at Holuhraun, Iceland, Derived from Airborne LiDAR and Photogrammetry

Roughness can be used to characterize the morphologies of a lava flow. It can be used to identify lava flow features, provide insight into eruption conditions, and link roughness pattern across a lava flow to emplacement conditions. In this study, we use both the topographic position index (TPI) and the one-dimensional Hurst exponent (H) to derive lava flow unit roughness on the 2014–2015 lava field at Holuhraun using both airborne LiDAR and photogrammetric datasets. The roughness assessment was acquired from four lava flow features: (1) spiny lava, (2) lava pond, (3) blocky surface, and (4) inflated channel. The TPI patterns on spiny lava and inflated channels show that the intermediate TPI values correspond to a small surficial slope indicating a flat and smooth surface. Lava pond is characterized by low to high TPI values and forms a wave-like pattern. Meanwhile, irregular transitions patterns from low to high TPI values indicate a rough surface that is found in blocky surface and flow margins. The surface roughness of these lava features falls within the H range of 0.30 ± 0.05 to 0.76 ± 0.04. The roughest surface is the blocky, and inflated lava flows appear to be the smoothest surface among these four lava units. In general, the Hurst exponent values in the 2014–2015 lava field at Holuhraun has a strong tendency in 0.5, both TPI and Hurst exponent successfully derive quantitative flow roughness.


Introduction
In the Earth Sciences, surface roughness is important for modeling natural phenomena and classifying features of interest [1,2]. Surface roughness refers to a topographic expression of surface profiles over various scales (i.e., centimeters, meters, kilometers) [1][2][3]. Quantitative approaches to estimate the roughness of natural materials are increasingly sought in modern geological research [1]. Statistical descriptors of surface morphology, or roughness, are found in many applications, including volcanology, especially for analyzing lava flows. Field observations have long been used in the study of surface roughness of lava flow [2][3][4][5][6]. These analyses are mostly based on in situ measurements, which require extended time in the field [2,[6][7][8]. In practice, a grid is laid out on the sample surface, and heights are measured manually or with a profiling instrument [2][3][4] and continuous Global Positioning System (GPS) [7]. Lava roughness reflects the morphology of lava flow, and some flows can be distinguished by roughness [2,7,9]. Thus, roughness can be used to identify lava flow features that reflect eruption conditions. Patterns of roughness across the lava flow can be tied to emplacement conditions such as rate of flow and viscosity [2,3,[10][11][12]. Roughness is typically determined from topography [2,3,13] or radar backscattering [4,14,15] (e.g., root mean square (RMS) height, correlation length, and autocorrelation function). Generally, the quantification of surface roughness is derived from analyzing height variations along profiles.
In the past 25 years, the study of lava flows roughness has embraced the use of remote sensing for acquiring surface geometry [11,16,17]. Airborne light detection and ranging (LiDAR) scanning and photogrammetry offers rapid three-dimensional (3D) data capture and have made datasets increasingly available to scientists [1]. Similar measurement techniques have been used to obtain ground-based measurements of surface roughness [7,16]. High-resolution topographic methods such as airborne LiDAR and photogrammetry generated DEMs (Digital Elevation Models) with meter scale resolution [8,18] allow for detailed roughness assessment. In this study, we evaluated two approaches to assess lava surface roughness based on LiDAR and photogrammetry datasets using: (1) the Topographic Position Index (TPI) [11,19]; and (2) the Hurst exponent (H) [2,7,8,16] on the most recent effusive eruption in Iceland, the 2014-2015 lava field at Holuhraun. We envisage future applications to assess geomorphic variation amongst different lava flows on Earth and other planets for which DEMs of high resolution and large areal extent are available.

The Surface Roughness of Lava Flows
Measurements of a lava flow surface roughness on Earth are used to describe changes in eruption conditions across a flow [20][21][22], surface processes that have occurred post-emplacement [2], to map flow units [3,8] and the relation between lava roughness and composition [2,9,23]. Effusive eruptions of basaltic magmas generally produce lava flows due to high magma temperatures and low viscosity. As lava starts to flow and cool, its surface may fold or break into blocks if the surface is steeply sloping due to the flow front stagnating from cooling [9,24]. This influences centimeter to decameter-scale roughness. If the emplacement surface is flat, the flow will advance slowly and take longer to cool, producing a smooth crust while intact, and a rougher surface if the crust breaks into blocks [9,11].
According to Kilburn [25], most basaltic lavas can be grouped according to surface roughness; (1) pahoehoe, where the surface is smooth and continuous (Figure 1a), (2) aa, where the surface is rough and fragmented (Figure 1b), and (3) blocky, where the surface is brecciated (Figure 1c). Furthermore, a transition surface morphology between aa and pahoehoe has been described as platty, slabby, and rubbly [26]. Thus, assessing lava flow textures can provide insight into the lava flow dynamics. Lower viscosities or shear strain results in smooth textures; rough textures are generally the result of higher viscosities, higher shear strain, or disruption of the cooled surface [2,3,7]. Identification of lava flows textures can further explain the geologic history of the eruption by confirming the styles, timing, and geographic extent of volcanic activity that occurred [27].
Volcanic eruptions can be understood through the variations in lava flow roughness at different scales. Spiny pahoehoe (Figure 1a) are typically smooth at meter scale and with spinose preserved on the surface of a pahoehoe flow and characterized by longitudinal grooves and ridges [28]. Ropy folds (Figure 1d) also provide an example of relating roughness to emplacement conditions [29,30]. Field observations of solidified pahoehoe surfaces and motion pictures of active flows suggest that these features may be interpreted as folds that develop in response to the shortening of the flow surface [30]. These features also indicate slowly moving low viscosity pahoehoe flows in basaltic eruptions [20,23,27,[30][31][32]. Though these features may appear flat and devoid of height changes at the meter and decameter scale, the texture is quite rough when observed at the centimeter scale. Spiny and ropy folds are both different from aa lava flows, which are composed of piles of jagged blocks and are rough at centimeter to decameter scale [2,3]. Table 1 shows lava surface features and the scales at which they are observed. Geosciences 2019, 9, x FOR PEER REVIEW 3 of 23  A variety of methods are used to quantify surface roughness of lava surfaces [2,3,35]. Two commonly used methods are the RMS of height and the H [2,7,16]. RMS represents the standard deviation of the height slope around the mean height [2]. RMS is a valuable method for vertical roughness, but it does not account for the horizontal patterns [11,36]. It has been used to study the roughness of lunar impact melts, and Martian lava flows [7]. To accurately reflect lava surface roughness, elevations in 360 degrees around a point should be considered, rather than only measuring topographic changes along one horizontal direction [11]. Therefore, we use TPI, and the one-dimensional H for determining the roughness of the 2014-2015 Holuhraun lava flows from both LiDAR DEM and photogrammetry DEM. The detailed methodology for these techniques will be explained in Section 4.2.  A variety of methods are used to quantify surface roughness of lava surfaces [2,3,35]. Two commonly used methods are the RMS of height and the H [2,7,16]. RMS represents the standard deviation of the height slope around the mean height [2]. RMS is a valuable method for vertical roughness, but it does not account for the horizontal patterns [11,36]. It has been used to study the roughness of lunar impact melts, and Martian lava flows [7]. To accurately reflect lava surface roughness, elevations in 360 degrees around a point should be considered, rather than only measuring topographic changes along one horizontal direction [11]. Therefore, we use TPI, and the one-dimensional H for determining the roughness of the 2014-2015 Holuhraun lava flows Geosciences 2020, 10, 125 4 of 23 from both LiDAR DEM and photogrammetry DEM. The detailed methodology for these techniques will be explained in Section 4.2.

Morphology of the 2014-2015 Lava Flow at Holuhraun
The six months long eruption at Holuhraun 2014-2015 was the largest effusive eruption in Iceland in 230 years with an estimated bulk lava volume of 1.44 km 3 [37][38][39]. The eruption is split into three phases based on the lava flow-field evolution [37]. The first phase was dominated by open channels, and during the second phase, lava emplacement was affected by the formation of a 1 km 2 lava pond about 1 km downstream of the vent [37]. This pond became the main distribution point for the lava during this phase, controlling the emplacement of lava flows. Near the end of the second phase, vertical stacking of lava lobes became more prevalent, and lava tubes developed within the channel system, resulting in the formation of inflation plateaus [37]. In the final phase, transport of lava through tubes continued, and inflation plateaus grew in extent, raising the original channel surface by 5-10 m above the surrounding lava. Over 19 km 2 of the flow field was resurfaced via surface breakouts from the closed paths during this period. Pedersen et al. [37] suggest that the topography of the lava field and surrounding made it possible to build an open channel system that was at minimum 5-10 m higher than the fluid lava. This system increased the static lava pressure, which was enough to lift the roof of the lava channels creating the inflation plateaus, allowing new lava to be transported to the distal ends of the lava field.
During the three phases of the 2014-2015 eruption, major changes in surface morphology of the lava occurred several times [37,40]. According to Pedersen et al. [37], shelly pahoehoe, slabby pahoehoe, rubbly pahoehoe, spiny pahoehoe, and aa were observed within the first week of the eruption. During the first phase and the second phase, aa was the dominant flow morphology of lava flows, and in the final phase, spiny pahoehoe was the main lava morphology [37]. This change from aa and pahoehoe morphology in the first and second phases to spiny pahoehoe in the final phase makes Holuhraun a paired lava flow-field. A paired lava flow-field is formed due to the decline of the effusion rate over the course of an eruption [41]. The Holuhraun lava flow-field was emplaced on a low-slope floodplain, and the chemical composition of the lava was uniform throughout the whole eruption [42,43]. This suggests that neither the topography nor the lava composition was the main factor for the observed changes in flow morphology [37]. The first transition of slabby pahoehoe to rubbly pahoehoe to aa occurred downstream of the vent, which is consistent with such changes in other lava producing eruptions and is explained by increased viscosity due to mixing, cooling, and gas loss during lava transport [34,37].
In this study, we examined roughness over lava flow surfaces in four locations depicted in Figure 2 that exhibit the known lava flow morphology at Holuhraun [7,37,44]. These lava flow features are: (1) lava pond, (2) spiny pahoehoe, (3) inflated channel, and (4) blocky surface. The lava pond formed during the first phase of the eruption [37]. Ropy folds preserved on the surface of a lava pond ( Figure 1d) indicate slowly moving, low viscosity pahoehoe flows, which develop in response to the shortening of the flow surface in basaltic eruptions [30]. At the Holuhraun lava flow-field, spiny lava is characterized by a network of interconnected lobes that form inflated sheet-like flow units with rough spinose surfaces [7]. The millimeter-scale spines on the surface of these flows resemble the texture of aa, but the flow surfaces are generally continuous and are not decomposed into clinker [7]. Spiny lava units are the dominant flow type along most of the flow margins, except near the NE margin of the flow [7]. The inflated channel formed towards the end of the second phase to the final phase [37]. This surface feature maintains the morphology of the flow channel but increases in thickness. Pahoehoe flows are typically identified as inflated lava flows, but aa type flows may also inflate under certain circumstances [45]. In this study, we use the term blocky surface to avoid confusion with rubbly lava. Blocky surface is composed of larger blocks than rubbly. Similar to rubbly lava this feature likely formed due to continued auto-brecciation of crustal slabs into blocks of material through mechanical collisions between the slabs during transport [26].At Holuhraun lava flow-field, Geosciences 2020, 10, 125 5 of 23 blocky surfaces are mostly found close to the vent. These lava flows were analyzed using the approach described in Section 4, and the results for the roughness analysis are presented in Section 5.

Datasets and Methods
To assess the surface roughness of lava flow features, we use high-DEM from two different sources. We focus on a subset of these data, obtained at four different facies on the lava flow-field. The datasets used in this study are described further in Table 2 and Section 4.1

Datasets and Methods
To assess the surface roughness of lava flow features, we use high-DEM from two different sources. We focus on a subset of these data, obtained at four different facies on the lava flow-field. The datasets used in this study are described further in Table 2 and Section 4.1

Airborne LiDAR and Photogrammetry
Airborne LiDAR data, collected and processed by the Natural Environment Research Council (NERC) was acquired on September 4, 2015. We processed the point clouds to obtain 1 m spatial resolution over the lava flow-field with a vertical resolution of 4-5 cm (depending on the flight line). This was done using ENVI LiDAR 5.3 software. Eight flight lines were flown over Holuhraun: seven of these are parallel and aligned with the long axis of the field, while the eighth is transverse and crosses all the others (Figure 3a). The LiDAR measurements, therefore, do not cover the entirety of the flow-field. Furthermore, small clouds and fumaroles obscured parts of the lava and created gaps in the data. We sought to determine the surface roughness of the 2014-2015 lava flows at Holuhraun. DEMs are produced upon processing of point-clouds for a variety of lava flow textures around the study area. Airborne LiDAR data, collected and processed by the Natural Environment Research Council (NERC) was acquired on September 4, 2015. We processed the point clouds to obtain 1 m spatial resolution over the lava flow-field with a vertical resolution of 4-5 cm (depending on the flight line). This was done using ENVI LiDAR 5.3 software. Eight flight lines were flown over Holuhraun: seven of these are parallel and aligned with the long axis of the field, while the eighth is transverse and crosses all the others (Figure 3a). The LiDAR measurements, therefore, do not cover the entirety of the flow-field. Furthermore, small clouds and fumaroles obscured parts of the lava and created gaps in the data. We sought to determine the surface roughness of the 2014-2015 lava flows at Holuhraun. DEMs are produced upon processing of point-clouds for a variety of lava flow textures around the study area.
We used photogrammetry-derived DEM provided by Loftmyndir ehf using data taken on August 30, 2015. Clouds obscured in some parts of the 2014-2015 Holuhraun lava flow-field ( Figure  2). This Photogrammetry-derived DEM has a spatial resolution of 5 m (Figure 3b).

Deriving Surface Roughness
The methodology adopted for the assessment of the roughness is given in a sequential manner in  We used photogrammetry-derived DEM provided by Loftmyndir ehf using data taken on August 30, 2015. Clouds obscured in some parts of the 2014-2015 Holuhraun lava flow-field ( Figure 2). This Photogrammetry-derived DEM has a spatial resolution of 5 m (Figure 3b).

Deriving Surface Roughness
The methodology adopted for the assessment of the roughness is given in a sequential manner in Figure 4. In this study, we use TPI and one-dimensional H to quantify roughness values in the 2014-2015 Holuhraun lava flows from both photogrammetry DEM and LiDAR DEM. Finally, we determined the roughness properties for four different facies on the lava flow-field from Figure 2. Airborne LiDAR data, collected and processed by the Natural Environment Research Council (NERC) was acquired on September 4, 2015. We processed the point clouds to obtain 1 m spatial resolution over the lava flow-field with a vertical resolution of 4-5 cm (depending on the flight line). This was done using ENVI LiDAR 5.3 software. Eight flight lines were flown over Holuhraun: seven of these are parallel and aligned with the long axis of the field, while the eighth is transverse and crosses all the others (Figure 3a). The LiDAR measurements, therefore, do not cover the entirety of the flow-field. Furthermore, small clouds and fumaroles obscured parts of the lava and created gaps in the data. We sought to determine the surface roughness of the 2014-2015 lava flows at Holuhraun. DEMs are produced upon processing of point-clouds for a variety of lava flow textures around the study area.
We used photogrammetry-derived DEM provided by Loftmyndir ehf using data taken on August 30, 2015. Clouds obscured in some parts of the 2014-2015 Holuhraun lava flow-field ( Figure  2). This Photogrammetry-derived DEM has a spatial resolution of 5 m (Figure 3b).

Deriving Surface Roughness
The methodology adopted for the assessment of the roughness is given in a sequential manner in Figure 4. In this study, we use TPI and one-dimensional H to quantify roughness values in the 2014-2015 Holuhraun lava flows from both photogrammetry DEM and LiDAR DEM. Finally, we determined the roughness properties for four different facies on the lava flow-field from Figure 2.

Slope
The estimation of slope from a regularly gridded DEM is a common procedure in terrain analysis [46,47]. The slope can be defined as a function of the gradients in the x and y direction at every point in a DEM: The key in the slope estimation is the computation of the perpendicular gradients f x and f y . We used moving 3 × 3 windows to derive local polynomial surface fit for the calculation [46,47]. This process was done using ArcGIS 10.6.

Topographic Position Index
In this study, topographic position index (TPI) [48] was used for deriving the roughness pattern of the lava flow-field. This method was originally created for use in ecology, geomorphology, and hydrology study [19,48] but has been recently used for assessing the topographic characteristics of lava flows [11]. TPI compares the elevation of each cell in a DEM to the mean elevation of a specified neighborhood around that cell. Mean elevation is subtracted from the elevation value at the center where C 0 is the elevation of the model point under evaluation, C is the mean elevation of a gridpoint in the neighborhood, σ is the standard deviation of elevation in the neighborhood. Positive TPI indicates that a cell is higher in elevation (or more steeply sloping) than the average of its neighbors up to a specified distance away, whereas a negative one indicates that a cell is lower than the average of surrounding elevations ( Figure 5) [19]. The cell neighborhood can be adapted to produce varying TPI values for different scales, thus changing the scale of roughness could affect the results [11,19,48].
In this study, we use a rectangular TPI with a 3 × 3 neighborhood size for both LiDAR DEM and photogrammetry DEM.

Slope
The estimation of slope from a regularly gridded DEM is a common procedure in terrain analysis [46,47]. The slope can be defined as a function of the gradients in the x and y direction at every point in a DEM: The key in the slope estimation is the computation of the perpendicular gradients and . We used moving 3 × 3 windows to derive local polynomial surface fit for the calculation [46,47]. This process was done using ArcGIS 10.6.

Topographic Position Index
In this study, topographic position index (TPI) [48] was used for deriving the roughness pattern of the lava flow-field. This method was originally created for use in ecology, geomorphology, and hydrology study [19,48] but has been recently used for assessing the topographic characteristics of lava flows [11]. TPI compares the elevation of each cell in a DEM to the mean elevation of a specified neighborhood around that cell. Mean elevation is subtracted from the elevation value at the center where is the elevation of the model point under evaluation, ̅ is the mean elevation of a gridpoint in the neighborhood, is the standard deviation of elevation in the neighborhood. Positive TPI indicates that a cell is higher in elevation (or more steeply sloping) than the average of its neighbors up to a specified distance away, whereas a negative one indicates that a cell is lower than the average of surrounding elevations ( Figure 5) [19]. The cell neighborhood can be adapted to produce varying TPI values for different scales, thus changing the scale of roughness could affect the results [11,19,48]. In this study, we use a rectangular TPI with a 3 × 3 neighborhood size for both LiDAR DEM and photogrammetry DEM.  [19]. The red lines indicate the neighborhood size.  [19]. The red lines indicate the neighborhood size.

Hurst Exponent
In this study, we acquired a one-dimensional TPI profile over various lava flow surfaces in four locations in Figure 2. The profiles were typically a few hundred meters to kilometers long. Roughness properties may have directional biases [2,7,8]. In this study, profiles were collected in perpendicular directions ( Figure 2). The Hurst exponent (H) is derived using rescaled range analysis (R/S) [15,49] for the TPI profiles. Equation (3) below shows the relation between R/S and H as where R is the difference between the maximum and minimum TPI detected in profiles, S is the standard deviation, and τ is the transect profiles. Hurst coefficient ranges from 0 to 1, where a higher H means a smooth or equally rough surface [2]. The mean Hurst exponent (H) was calculated for each lava unit.

Results
We extracted the slope, TPI, and Hurst exponent from the lava flow features in Figure 2 (lava pond, spiny pahoehoe, inflated channel, and blocky surface). Each lava feature displays TPI and H variations, although the entire range of H variation is not displayed at any individual flow. The roughness analysis results are shown as maps and profiles for each lava flow features. Based on the analysis, the surface feature that corresponds to a higher H and intermediate TPI pattern reflects a smoother surface than the lower H and irregular TPI patterns.

Lava Pond Roughness
For the photogrammetry DEM (Figure 6a), TPI and slope values at the lava pond ranged from −2.8 m to 2.8 m and from 0.4 • to 12.5 • (Figure 6b,c), respectively. For the LiDAR DEM (Figure 6d), TPI and slope values ranged from -3.2 m to 5.4 m and from 0.5 • to 30 • (Figure 6e,f), respectively. Both the LiDAR and photogrammetry derived TPI images have wave-like transitions patterns (Figure 6c,f). These patterns are also apparent from the slope (Figure 6b,e), indicating a relatively rough surface due to the ropy fold on the lava pond. Folds form where the velocity and viscosity of the lava surface decreases rapidly with depth [30]. This causes surface folding and irregular waviness in the surface which preserve on wave-like pattern on TPI maps. In this study, we acquired a one-dimensional TPI profile over various lava flow surfaces in four locations in Figure 2. The profiles were typically a few hundred meters to kilometers long. Roughness properties may have directional biases [2,7,8]. In this study, profiles were collected in perpendicular directions ( Figure 2). The Hurst exponent (H) is derived using rescaled range analysis (R/S) [15,49] for the TPI profiles. Equation (3)

below shows the relation between R/S and H as
where R is the difference between the maximum and minimum TPI detected in profiles, S is the standard deviation, and τ is the transect profiles. Hurst coefficient ranges from 0 to 1, where a higher H means a smooth or equally rough surface [2]. The mean Hurst exponent (H ) was calculated for each lava unit.

Results
We extracted the slope, TPI, and Hurst exponent from the lava flow features in Figure 2 (lava pond, spiny pahoehoe, inflated channel, and blocky surface). Each lava feature displays TPI and H variations, although the entire range of H variation is not displayed at any individual flow. The roughness analysis results are shown as maps and profiles for each lava flow features. Based on the analysis, the surface feature that corresponds to a higher H and intermediate TPI pattern reflects a smoother surface than the lower H and irregular TPI patterns.

Lava Pond Roughness
For the photogrammetry DEM (Figure 6a), TPI and slope values at the lava pond ranged from −2.8 m to 2.8 m and from 0.4° to 12.5° (Figure 6b,c), respectively. For the LiDAR DEM (Figure 6d), TPI and slope values ranged from -3.2 m to 5.4 m and from 0.5° to 30° (Figure 6e,f), respectively. Both the LiDAR and photogrammetry derived TPI images have wave-like transitions patterns (Figure 6c,f). These patterns are also apparent from the slope (Figure 6b,e), indicating a relatively rough surface due to the ropy fold on the lava pond. Folds form where the velocity and viscosity of the lava surface decreases rapidly with depth [30]. This causes surface folding and irregular waviness in the surface which preserve on wave-like pattern on TPI maps.

Spiny Lava Roughness
The photogrammetry DEM (Figure 7a) had TPI and slope values at the spiny lava ranging from -1.8 m to 2.2 m and from 0° to 38° (Figure 7b,c), respectively. On the other hand, the LiDAR DEM (Figure 7d), had TPI and slope values ranging from -2.1 m to 1.8 m and from 0° to 36° ( Figure  6e,f). The LiDAR DEM does not cover the entirety of the spiny lava due to gaps in the image. The low TPI values correspond to a low slope and indicate a flat and smooth surface that has formed from inflated spiny lava; however, we cannot differentiate the spinose pattern in this meter scale. The irregular transitions from low to high TPI and slope pattern around the flow margin might indicate a rough surface resulting from a breakout from the inflated spiny lava. These features will be explained in more detail in Section 6.2.

Spiny Lava Roughness
The photogrammetry DEM (Figure 7a) had TPI and slope values at the spiny lava ranging from -1.8 m to 2.2 m and from 0 • to 38 • (Figure 7b,c), respectively. On the other hand, the LiDAR DEM (Figure 7d), had TPI and slope values ranging from -2.1 m to 1.8 m and from 0 • to 36 • (Figure 6e,f). The LiDAR DEM does not cover the entirety of the spiny lava due to gaps in the image. The low TPI values correspond to a low slope and indicate a flat and smooth surface that has formed from inflated spiny lava; however, we cannot differentiate the spinose pattern in this meter scale. The irregular transitions from low to high TPI and slope pattern around the flow margin might indicate a rough surface resulting from a breakout from the inflated spiny lava. These features will be explained in more detail in Section 6.2.

Spiny Lava Roughness
The photogrammetry DEM (Figure 7a) had TPI and slope values at the spiny lava ranging from -1.8 m to 2.2 m and from 0° to 38° (Figure 7b,c), respectively. On the other hand, the LiDAR DEM (Figure 7d), had TPI and slope values ranging from -2.1 m to 1.8 m and from 0° to 36° ( Figure  6e,f). The LiDAR DEM does not cover the entirety of the spiny lava due to gaps in the image. The low TPI values correspond to a low slope and indicate a flat and smooth surface that has formed from inflated spiny lava; however, we cannot differentiate the spinose pattern in this meter scale. The irregular transitions from low to high TPI and slope pattern around the flow margin might indicate a rough surface resulting from a breakout from the inflated spiny lava. These features will be explained in more detail in Section 6.2.

Inflated Channel Roughness
The photogrammetry DEM at the inflated channel (Figure 8a) yielded TPI and slope ranging from -3 m to 3.2 m and from 0.1° to 25° (Figure 8b,c), respectively. The LiDAR DEM (Figure 8d) yielded TPI and slope values ranging from -11 m to 4.8 m and from 0° to 33° (Figure 8e,f). Comparable to what we have found in spiny lava, low TPI patterns that correspond to a low slope were assigned to the flat surface of the inflated channel. The lowest TPI values correspond to a high slope indicating inflation pits and cracks. Inflation pits form where a portion of the flow is not inflating. The highest TPI values correspond to a high slope relating to boulders and grooves within the lava. The irregular transitions pattern from a low to a high slope and TPI values around the inflation margin are similar to what we have found on spiny lava. This indicates this feature is a rough surface that results from a breakout from the inflated channel (see Section 6.2 for more detail).

Inflated Channel Roughness
The photogrammetry DEM at the inflated channel (Figure 8a) yielded TPI and slope ranging from -3 m to 3.2 m and from 0.1 • to 25 • (Figure 8b,c), respectively. The LiDAR DEM (Figure 8d) yielded TPI and slope values ranging from -11 m to 4.8 m and from 0 • to 33 • (Figure 8e,f). Comparable to what we have found in spiny lava, low TPI patterns that correspond to a low slope were assigned to the flat surface of the inflated channel. The lowest TPI values correspond to a high slope indicating inflation pits and cracks. Inflation pits form where a portion of the flow is not inflating. The highest TPI values correspond to a high slope relating to boulders and grooves within the lava. The irregular transitions pattern from a low to a high slope and TPI values around the inflation margin are similar to what we have found on spiny lava. This indicates this feature is a rough surface that results from a breakout from the inflated channel (see Section 6.2 for more detail).

Inflated Channel Roughness
The photogrammetry DEM at the inflated channel ( Figure 8a) yielded TPI and slope ranging from -3 m to 3.2 m and from 0.1° to 25° (Figure 8b,c), respectively. The LiDAR DEM (Figure 8d) yielded TPI and slope values ranging from -11 m to 4.8 m and from 0° to 33° (Figure 8e,f). Comparable to what we have found in spiny lava, low TPI patterns that correspond to a low slope were assigned to the flat surface of the inflated channel. The lowest TPI values correspond to a high slope indicating inflation pits and cracks. Inflation pits form where a portion of the flow is not inflating. The highest TPI values correspond to a high slope relating to boulders and grooves within the lava. The irregular transitions pattern from a low to a high slope and TPI values around the inflation margin are similar to what we have found on spiny lava. This indicates this feature is a rough surface that results from a breakout from the inflated channel (see Section 6.2 for more detail).

Block Surface Roughness
The photogrammetry DEM (Figure 9a) at the blocky surface units had TPI and slope ranging from -2.1 m to 2 m and from 0.4° to 15°, respectively (Figure 9a,b). The LiDAR DEM (Figure 9d) had TPI, and slope values ranging from −2.3 m to 2.3 m and from 0° to 23°, respectively (Figure 9e,f). In the LiDAR, blocky surface, typically, had irregular TPI patterns. These patterns indicate that the surface is rough which related to the decimeter to meter scale blocks of fragmented pahoehoe-like crust. However, in the photogrammetry DEM, blocky surface appears to have smoother patterns corresponding to high slope since the roughness of the blocky surface is not distinct in lower resolution.

Block Surface Roughness
The photogrammetry DEM (Figure 9a) at the blocky surface units had TPI and slope ranging from -2.1 m to 2 m and from 0.4 • to 15 • , respectively (Figure 9a In the LiDAR, blocky surface, typically, had irregular TPI patterns. These patterns indicate that the surface is rough which related to the decimeter to meter scale blocks of fragmented pahoehoe-like crust. However, in the photogrammetry DEM, blocky surface appears to have smoother patterns corresponding to high slope since the roughness of the blocky surface is not distinct in lower resolution.

Block Surface Roughness
The photogrammetry DEM (Figure 9a) at the blocky surface units had TPI and slope ranging from -2.1 m to 2 m and from 0.4° to 15°, respectively (Figure 9a,b). The LiDAR DEM (Figure 9d) had TPI, and slope values ranging from −2.3 m to 2.3 m and from 0° to 23°, respectively (Figure 9e,f). In the LiDAR, blocky surface, typically, had irregular TPI patterns. These patterns indicate that the surface is rough which related to the decimeter to meter scale blocks of fragmented pahoehoe-like crust. However, in the photogrammetry DEM, blocky surface appears to have smoother patterns corresponding to high slope since the roughness of the blocky surface is not distinct in lower resolution.

One dimensional TPI profiles
The one-dimensional TPI profiles over the four lava flow units are shown in Figure 10a-d. Based on H reported in Table 3, the surface roughness of lava features falls within the range from 0.3 ± 0.05 to 0.76 ± 0.04. Typically, LiDAR DEM yields lower Hurst exponent values than photogrammetry DEM. This was also found for blocky surface, i.e., photogrammetry DEM had higher Hurst exponent than LiDAR DEM. LiDAR DEM has greater pixel resolution than photogrammetry DEM which results in more detailed TPI profiles. This issue will be addressed in Section 6.1. Based on the Mean Hurst exponent (H ), the roughest surface is the blocky surface with H ' around 0.52 ± 0.04 and the inflated lava field appears to be the smoothest surface of these four lava units with H around 0.61 ± 0.06. These results are comparable with the results from the TPI map pattern that showed that inflated lava appears to be the smoothest surface, and blocky is the roughest surface. In general, the Hurst exponent values have a strong tendency to be close to 0.5. This was also suggested by an early study by Shepard et al. [2] for geological surface roughness.

One Dimensional TPI Profiles
The one-dimensional TPI profiles over the four lava flow units are shown in Figure 10a-d. Based on H reported in Table 3, the surface roughness of lava features falls within the range from 0.3 ± 0.05 to 0.76 ± 0.04. Typically, LiDAR DEM yields lower Hurst exponent values than photogrammetry DEM. This was also found for blocky surface, i.e., photogrammetry DEM had higher Hurst exponent than LiDAR DEM. LiDAR DEM has greater pixel resolution than photogrammetry DEM which results in more detailed TPI profiles. This issue will be addressed in Section 6.1. Based on the Mean Hurst exponent (H), the roughest surface is the blocky surface with H' around 0.52 ± 0.04 and the inflated lava field appears to be the smoothest surface of these four lava units with H around 0.61 ± 0.06. These results are comparable with the results from the TPI map pattern that showed that inflated lava appears to be the smoothest surface, and blocky is the roughest surface. In general, the Hurst exponent values have a strong tendency to be close to 0.5. This was also suggested by an early study by Shepard et al. [2] for geological surface roughness.

Discussion
In this work, we examined lava flows roughness in the 2014-2015 lava field at Holuhraun using TPI and one-dimensional H profiles. Both TPI and H successfully derived lava roughness on selected lava units at Holuhraun. However, there are still some issues. How does the scale affect the roughness results for both TPI and Hurst exponent? Is the 'smooth' surface still smooth for a different scale? How do we connect the roughness to the emplacement style? What can we improve to get a better roughness assessment? In this section, we address several issues related to (1) TPI neighborhood size and profile length, (2) Eruption condition and link roughness pattern with the emplacement style, and (3) alternative datasets and methods for deriving roughness.

TPI Neighborhood Size and Profile Length
TPI is naturally very scale dependent [11,19,48,50]. We should consider what scale is the most relevant for the study being analyzed [19]. The patterns produced by TPI vary on the scale that we use to analyze. Figure 11 shows an illustration of similar topography using different neighborhood sizes in TPI [19,51]. We test different TPI neighborhood sizes on LiDAR DEM spiny lava, as shown in Figure 12a-d. The 1 × 1 neighborhood size (Figure 12a) was not as clear as the larger scales to differentiate flow margins but appears to emphasize small features on the lava flows, which does not show in larger neighborhoods. James [11] recommends building a catalog of features presented at each scale. This could be useful in assessing the topographic characteristics of a lava flow surface. relevant for the study being analyzed [19]. The patterns produced by TPI vary on the scale that we use to analyze. Figure 11 shows an illustration of similar topography using different neighborhood sizes in TPI [19,51]. We test different TPI neighborhood sizes on LiDAR DEM spiny lava, as shown in Figure 12a-d. The 1 × 1 neighborhood size (Figure 12a) was not as clear as the larger scales to differentiate flow margins but appears to emphasize small features on the lava flows, which does not show in larger neighborhoods. James [11] recommends building a catalog of features presented at each scale. This could be useful in assessing the topographic characteristics of a lava flow surface. Figure 11. The illustration of the effects of the neighborhood size on TPI value adapted and modified from Jennes [19]. The red line depicts the neighborhood size. In topography (A), the neighborhood size is small enough that the point is at about the same elevation as the entire analysis region, so the TPI value is approximately 0, which means it is a flat surface. In topography (B), the neighborhood size is large enough to encompass the entire small hill, and the point is consequently much higher than its surroundings and has a correspondingly high TPI value resulting in that the point is detected as a ridge [19]. In topography (C), the neighborhood includes the hills on either side of the valley, and therefore the point is lower than its neighbors and has a negative TPI value [11,19].
Geosciences 2019, 9, x FOR PEER REVIEW 15 of 23 Figure 11. The illustration of the effects of the neighborhood size on TPI value adapted and modified from Jennes [19]. The red line depicts the neighborhood size. In topography A, the neighborhood size is small enough that the point is at about the same elevation as the entire analysis region, so the TPI value is approximately 0, which means it is a flat surface. In topography B, the neighborhood size is large enough to encompass the entire small hill, and the point is consequently much higher than its surroundings and has a correspondingly high TPI value resulting in that the point is detected as a ridge [19]. In topography C, the neighborhood includes the hills on either side of the valley, and therefore the point is lower than its neighbors and has a negative TPI value [11,19]. We consider that there are at least two factors that affect the H values: (1) pixel size and (2) the profile length. Figure 13a shows TPI profiles representing blocky surface vertical profile (BLV1) from LiDAR DEM and photogrammetry DEM. The LiDAR DEM profile is rougher than the photogrammetry DEM since the spatial spacing (pixel size) in the LiDAR DEM is denser than photogrammetry DEM because the resolution is higher (1 m vs. 5 m). This spatial spacing issue corresponds to the TPI neighborhood size in Figure 11, where LiDAR DEM captures more detailed TPI value than photogrammetry DEM which affects the Hurst exponent derived roughness. Second, we also consider the profile length as a factor that affects the Hurst exponent values. In Figure 13b, we subset the first 100 m BLV1 profiles, and the Hurst exponent values for these subset profiles were increased for both topographies. For LiDAR DEM the Hurst exponent values are 0.55 ± 0.05 compared to 0.30 ± 0.05 for full profiles, and for photogrammetry DEM the Hurst exponent values are 0.48 ± 0.01 compared to 0.40 ± 0.05 for full profiles. These results not necessarily conclude that the shorter profiles will increase the Hurst exponent. Extended profiles also can increase the Hurst exponent since the profiles could be equally rough, which can represent a 'smooth' surface [2]. In many cases, it depends on the morphology that we are studying [2,3], and we should know the morphology. The We consider that there are at least two factors that affect the H values: (1) pixel size and (2) the profile length. Figure 13a shows TPI profiles representing blocky surface vertical profile (BLV1) from LiDAR DEM and photogrammetry DEM. The LiDAR DEM profile is rougher than the photogrammetry DEM since the spatial spacing (pixel size) in the LiDAR DEM is denser than photogrammetry DEM because the resolution is higher (1 m vs. 5 m). This spatial spacing issue corresponds to the TPI neighborhood size in Figure 11, where LiDAR DEM captures more detailed TPI value than photogrammetry DEM which affects the Hurst exponent derived roughness. Second, we also consider the profile length as a factor that affects the Hurst exponent values. In Figure 13b, we subset the first 100 m BLV1 profiles, and the Hurst exponent values for these subset profiles were increased for both topographies. For LiDAR DEM the Hurst exponent values are 0.55 ± 0.05 compared to 0.30 ± 0.05 for full profiles, and for photogrammetry DEM the Hurst exponent values are 0.48 ± 0.01 compared to 0.40 ± 0.05 for full profiles. These results not necessarily conclude that the shorter profiles will increase the Hurst exponent. Extended profiles also can increase the Hurst exponent since the profiles could be equally rough, which can represent a 'smooth' surface [2]. In many cases, it depends on the morphology that we are studying [2,3], and we should know the morphology. The other thing we found is that the direction of profiles could bias the Hurst exponent [7,52]. We recommend in a future study to build series of profiles that are rotated by some number of degrees to capture a wider range of directions around the surface area of interest to derive precise roughness. Some studies are currently developing three-dimensional roughness estimates based on a 3D Gaussian filter applied to DEM [53].
Geosciences 2019, 9, x FOR PEER REVIEW 16 of 23 recommend in a future study to build series of profiles that are rotated by some number of degrees to capture a wider range of directions around the surface area of interest to derive precise roughness. Some studies are currently developing three-dimensional roughness estimates based on a 3D Gaussian filter applied to DEM [53].
(a) (b) Figure 13. TPI profiles representing the first vertical profile of the blocky surface (BLV1) from the LiDAR DEM (green) and the photogrammetry DEM (red). (a) full profiles; (b) a zoom-in of the first 100 meters of (a). The profiles had been offset from each other for clarity.

Eruption Condition and Link Quantitative Roughness with the Emplacement Style
Quantitative roughness measurement of lava flow features provides insight into the relative emplacement styles. In particular, higher H and intermediate TPI patterns suggest a smoother surface, which reflects relatively lower viscosity lavas than for the lower H and irregular TPI pattern. However, roughness variations do not directly distinguish between lava viscosity. To quantify precise viscosity, more factors need to consider rather than single roughness. The difference in viscosity may be the result of transport within an insulated distributary system that limits heat loss. If the differences between the flow types are predominantly textural, local supply may be a controlling factor [5].
The two primary roughness identified in the, a smooth surface possibly inflated flow type and a rough surface possibly breakout, blocky surface and ropy fold data. Section 5.2 and 5.3 shows that smoother surface types exhibit morphologic properties attributed to flow inflation. Along margins, irregular transitions TPI pattern supports the interpretation of the occurrence of breakout (Figure 7c-Figure 13. TPI profiles representing the first vertical profile of the blocky surface (BLV1) from the LiDAR DEM (green) and the photogrammetry DEM (red). (a) full profiles; (b) a zoom-in of the first 100 meters of (a). The profiles had been offset from each other for clarity.

Eruption Condition and Link Quantitative Roughness with the Emplacement Style
Quantitative roughness measurement of lava flow features provides insight into the relative emplacement styles. In particular, higher H and intermediate TPI patterns suggest a smoother surface, which reflects relatively lower viscosity lavas than for the lower H and irregular TPI pattern. However, roughness variations do not directly distinguish between lava viscosity. To quantify precise viscosity, more factors need to consider rather than single roughness. The difference in viscosity may be the result of transport within an insulated distributary system that limits heat loss. If the differences between the flow types are predominantly textural, local supply may be a controlling factor [5].
The two primary roughness identified in the, a smooth surface possibly inflated flow type and a rough surface possibly breakout, blocky surface and ropy fold data. Sections 5.2 and 5.3 shows that smoother surface types exhibit morphologic properties attributed to flow inflation. Along margins, irregular transitions TPI pattern supports the interpretation of the occurrence of breakout (Figures 7c-f and 8c-f). As displayed clearly in Figure 14a-b, the rough surface in lava flow margins from the three-dimension LiDAR surface model. These rough surfaces result from a breakout from the inflated surface from both spiny lava (Figure 14a) and inflated channel (Figure 14b). It is reported by Pedersen et al. [37] that a series of breakouts occurred during the third phase of the eruption, mostly in the flow margin and also in inflated channels. This breakout is inferred from the false color on the Earth Observing 1 (EO-1) Advanced Land Imager (ALI) satellite during the eruption (Figure 15). These breakouts occurred around 16 January 2015, a breakout from spiny lava and inflated channel. A breakout is a lobe originating from the liquid interior of active lava. It may take place through a crack at the front or the side of the flow margin when lava is inflated [54]. Injection of fresh basaltic magma into an established flow produces new breakouts of lava and promotes lifting of the upper crust (inflation). These breakout phenomena illustrated in Figure 16. inflated surface from both spiny lava (Figure 14a) and inflated channel (Figure 14b). It is reported by Pedersen et al. [37] that a series of breakouts occurred during the third phase of the eruption, mostly in the flow margin and also in inflated channels. This breakout is inferred from the false color on the Earth Observing 1 (EO-1) Advanced Land Imager (ALI) satellite during the eruption ( Figure 15). These breakouts occurred around 16 January 2015, a breakout from spiny lava and inflated channel. A breakout is a lobe originating from the liquid interior of active lava. It may take place through a crack at the front or the side of the flow margin when lava is inflated [54]. Injection of fresh basaltic magma into an established flow produces new breakouts of lava and promotes lifting of the upper crust (inflation). These breakout phenomena illustrated in Figure 16.

Alternative Datasets and Methods for Deriving Roughness
Recently, a number of and methods are used to quantify surface roughness. Neish et al. [7] quantified the surface roughness from synthetic aperture radar (SAR) using the circular polarization ratio (CPR). This technique is defined as the radar backscatter ratio in the same polarization that was transmitted (SC) to the opposite polarization (OC), where the CPR is SC divided by OC [7]. Rough surfaces tend to have approximately equal OC and SC returns, with CPR values approaching one. Meanwhile, the flat surfaces tend to have high OC returns and low CPR values [7]. We can also apply Hurst exponent to the radar backscattering profiles in order to derive roughness [15]. Neish et al. [7] also show that linking radar backscatter to the high-resolution topographic profiles is essential to understanding the structure of various geologic units. Another technique we could consider in describing roughness is optical images (multispectral-hyperspectral). Several studies characterize roughness from multispectral and hyperspectral based on the spectral reflectance [27,31,55,56]. The main assumption is that rough surface has lower reflectance than a smooth surface. The integration of these techniques with field measurement could improve surface roughness estimation of lava flow.

Conclusions
In this study, both TPI and one-dimensional H successfully derive quantitative flow roughness. The roughness assessment was acquired from four lava flow features: (1) spiny lava, (2) lava pond, (3) blocky surface, and (4) inflated channel. TPI patterns on spiny lava and inflated channels show that the intermediate TPI patterns with low slope indicate flat and smooth surface. Lava pond has transitions pattern from low to high TPI values forming a wave-like pattern. Meanwhile, irregular transitions from low to high TPI and the slope patterns indicate a rough surface that is found in blocky surface and flow margins, indicating lava breakouts. Quantitative measures of surface roughness of lava features fall within the H ranges from 0.30 ± 0.05 to 0.76 ± 0.04. The roughest surface is the blocky surface, and the inflated lava flow appears to be the smoothest surface among these four lava units. In general, the Hurst exponent values in the 2014-2015 lava field at Holuhraun has a strong tendency to be close to 0.5, which has good agreement with an earlier study for geological surface roughness. Neighborhood size is a critical component for TPI to quantify the roughness. Small neighborhoods capture small and local features and valleys, while large neighborhoods capture larger-scale features. We consider that there are at least two factors that affect the Hurst exponent values: (1) pixel size and (2) the profile length. We recommend in a future study to build a series of profiles that are rotated by some number of degrees to capture a wider range of directions, in order to derive precise roughness. The integration of multimodal remote sensing datasets and field measurement can improve the estimation of surface roughness of lava flow.