# Mapping Pyroclastic Flow Inundation Using Radar and Optical Satellite Images and Lahar Modeling

Authors: Chang-Wook Lee; Zhong Lu; Jin-Woo Kim; Sung-Jae Park
Journal: Journal of Sensors
 (2018)

---

## Abstract

Sinabung volcano, located above the Sumatra subduction of the Indo-Australian plate under the Eurasian plate, became active in 2010 after about 400 years of quiescence. We use ALOS/PALSAR interferometric synthetic aperture radar (InSAR) images to measure surface deformation from February 2007 to January 2011. We model the observed preeruption inflation and coeruption deflation using Mogi and prolate spheroid sources to infer volume changes of the magma chamber. We interpret that the inflation was due to magma accumulation in a shallow reservoir beneath Mount Sinabung and attribute the deflation due to magma withdrawal from the shallow reservoir during the eruption as well as thermoelastic compaction of erupted material. The pyroclastic flow extent during the eruption is then derived from the LAHARZ model based on the coeruption volume from InSAR modeling and compared to that derived from the Landsat 7 Enhanced Thematic Mapper Plus (ETM+) image. The pyroclastic flow inundation extents between the two different methods agree at about 86%, suggesting the capability of mapping pyroclastic flow inundation by combing radar and optical imagery as well as flow modeling.


---

## Body

## 1. Introduction

Volcanic pyroclastic flows consist of a large amount of destructive mass of very hot ash, lava fragments, and gases ejected explosively from a volcano and typically flow downslope at great speed. Volcanic pyroclastic flows pose great threats to people and property downstream from volcanoes. In order to accurately estimate the extents that are likely to be impacted by pyroclastic flows, computer models, such as LAHARZ [1], have been developed on the basis of statistical analyses of flow-path geometries and mathematical analyses of flow physics. LAHARZ [1] has been applied at many volcanoes worldwide.Mount Sinabung is a stratovolcano in the northwestern part of the Sumatra Island, Indonesia (Figure1). Formed during the Pleistocene to Holocene, the volcano consists of andesitic and dacitic lavas. The Sumatra Island, located on the boundary between the Australian Plate to the south and the Eurasian Plate to the north, was formed by the melting of subducted rocks under enormous pressure and its subsequent eruptions in volcanoes. Sumatra and Java, Indonesia, include a long line of volcanoes along the boundary between the two plates and have produced a large number of earthquakes and volcanic eruptions. Located within 50 km southeast of Sinabung, Sumatra Island’s Toba volcano had the largest eruption with a volcanic explosivity index (VEI) of 8 in the geological record during last 25 million years [2].Figure 1
Location of Mount Sinabung in Indonesia above the Sumatra subduction zone.Mount Sinabung has a 2460 m high conical structure with four overlapping craters about 60 to 300 m in diameter from north to south, respectively, and is mostly covered with dense vegetation exclusively at the top of the volcano (Figure2). Sinabung volcano had its last historical activity about 400 years ago, and the 2010 eruption started on August 29 and lasted to the middle of September 2010. After the first reawakened volcanic eruption in 2010, Mount Sinabung erupted again in 2013, 2014, and 2016.Figure 2
Landsat 7 ETM+ imagery at Mount Sinabung on 19 May 2003.SAR interferometry (InSAR) is an efficient method for measuring surface deformation over large areas [3–7], and it has been applied worldwide to study volcanoes [8–13]. Chaussard and Amelung [10] processed 13 L-band ALOS/PALSAR images acquired before 2010 and showed progressive inflation at a rate of 2.2 cm/yr using small-baseline subset (SBAS) InSAR technique [14–16]. They also measured the surface displacement before the eruption and during the eruption and estimated the depth of magma storage (0.9 km) using volcanic source modeling (i.e., Mogi and Yang) [17–21]. However, these previous InSAR studies did not calculate the substantive damages of the volcanic activities but merely observed the volume change of the magma chamber according to the arithmetical deformation of the surface.In this study, we have carried out a study of pyroclastic flow mapping using satellite imagery from both radar and optical sensors as well as the LAHARZ program (Figure3). We improve the preliminary results [22–24] by laying out details on how to derive the pyroclastic flow extent using maximum likelihood classification method based on Landsat 7 images where artifacts due to scan-line corrector have been rectified. We generate the coeruption InSAR image by stacking several interferograms in order to improve the signal-to-noise ratio (SNR) of the coeruption deformation map, based on which we derive the extrusion volume associated with the 2010 eruption at Sinabung from geophysical models. The goal is to generate a high-SNR deformation map for an accurate estimate of the eruption volume and improve the earlier results of time-series deformation measurements where noise was apparent in some individual InSAR images [22–24]. We derive the volume change associated with Sinabung’s volcanic eruption in 2010 using mechanical modeling, and the calculated volume change of the magma chamber is then applied to a pyroclastic flow modeling algorithm based on LAHARZ [1, 25–27] to obtain the hazard zone of pyroclastic flow deposit. We describe in detail how we model the pyroclastic inundation map using the LAHARZ model and discuss quantitatively the comparison of pyroclastic flow simulation maps and pyroclastic flow map derived from the Landsat 7 ETM+ image of Sinabung volcano (Table 1). Finally, we highlight that the integration of satellite image analysis and pyroclastic flow simulation using the LAHARZ program can accomplish volcano monitoring and pyroclastic flow inundation hazard mapping of active volcanoes and discuss the limitation of the proposed method.Figure 3
Flowchart of the radar and optical data processing.Table 1
Characteristics of ALOS/PALSAR and Landsat data used in this study.
(a)
Number
Mission
Orbit number
Date (YYYYMMDD)
Perpendicular baseline (m)
1
ALOS/PALSAR
00702
20070220
0
2
ALOS/PALSAR
00707
20070708
181
3
ALOS/PALSAR
00708
20070823
205
4
ALOS/PALSAR
00801
20080108
15
5
ALOS/PALSAR
00802
20080223
194
6
ALOS/PALSAR
00804
20080409
305
7
ALOS/PALSAR
00805
20080525
−45
8
ALOS/PALSAR
00810
20081010
−124
9
ALOS/PALSAR
00901
20090110
101
10
ALOS/PALSAR
00902
20090225
−227
11
ALOS/PALSAR
00907
20090713
449
12
ALOS/PALSAR
00908
20090828
546
13
ALOS/PALSAR
00911
20091128
147
14
ALOS/PALSAR
01001
20100113
−3
15
ALOS/PALSAR
01002
20100228
189
16
ALOS/PALSAR
01007
20100716
46
17
ALOS/PALSAR
01008
20100831
202
18
ALOS/PALSAR
01010
20101016
154
19
ALOS/PALSAR
01012
20101201
−165
20
ALOS/PALSAR
01101
20110116
31
(b)
Number
Mission
Path/Track
Date
Time (UTC)
1
Landsat/ETM+
129/058
20120730
03:35:52.8


## 2. Processing of Optical Image

A pyroclastic flow inundation map can be generated by field surveys from the crater rim to the furthest extent of the pyroclastic flow after a volcanic eruption. However, a field survey in an active volcano can be dangerous due to the exposure to hazardous gases and abrupt activities. In contrast, remote sensing techniques are a useful tool for generating pyroclastic flow deposit maps, providing a safe, cost effective alternative to field mapping.In this study, we aim to generate a pyroclastic flow deposit map using a Landsat 7 ETM+ image acquired on 30 Jul. 2012 after the 2010 Sinabung eruption (Table1). This image represents the land cover after the 2010 eruptions with an absence of atmospheric effects (Figure 4(a)). However, the Landsat 7 ETM+ data suffer from an artifact due to the failure of the scan-line corrector (SLC) in the ETM+ instrument. The Landsat 7 data can be restored by a gap-filling and filtering technique described by Lee et al. [28, 29] (Figure 4(b)). Finally, we use the gap-filled Landsat 7 ETM+ data (Figure 4(b)) to derive land-cover types through supervised maximum likelihood classification (Figure 4(c)). We establish six classes, including farm land (orange), forest (green), city (gray), shadow (dark gray), water (blue), and pyroclastic flow deposit area (red). Based on the classification maps, the pyroclastic flow deposits from the 2010 eruption were mainly contained over the southern and eastern upper flanks of the volcano (red in Figure 4(c)). The pyroclastic flow inundation area after the 2010 eruption of Sinabung volcano is estimated as 0.872 km2 (red in Figure 4(d)).Figure 4
(a) False color composite (RGB = 742) image of Landsat 7 ETM+ acquired on 30 Jul. 2012. (b) The gap-filled Landsat 7 ETM+ image. (c) Supervised classification map derived from (b). (d) Pyroclastic flow inundation after the 2010 eruption at Sinabung volcano.
(a)
(b)
(c)
(d)


## 3. Processing of Radar Images

InSAR processing is carried out to characterize the surface deformation of the Sinabung volcano. Using 20 ALOS/PALSAR images from 20 Feb. 2007 to 16 Jan. 2011 (Table1), 40 interferograms are generated. The best coherence was constrained on the volcano summit. Outside the summit, coherence was maintained only at many scattered patches of rocky and less-vegetated areas.The trend of surface deformation shows an uplift pattern until the end of August 2010 but a deflation pattern from the eruption event to the end of the dataset, corresponding to the changes in magma volume inside the reservoir. Therefore, we divide our results into two different periods of inflation and deflation separated by the onset of the 2010 eruption. The time-series SBAS technique is applied to measure the surface deformation of the Sinabung volcano at each period [16]. Using this method, we can minimize artificial atmospheric effects [30, 31] and accurately measure surface deformation over a specific time period. Figure 5(a) displays the stacked image showing inflation of magma volume change until the 2010 eruption, and Figure 5(b) represents the stacked image that shows the deflation pattern of surface deformation during the 2010 eruption.Figure 5
Stacked unwrapped interferograms for preeruption (a) and coeruption (b) of Sinabung volcano (modified from [16]).
(a)
(b)The two stacked images of surface deformation (Figures6(a) and 6(b)) have opposite fringe tendency, with a maximum 1.2 cm/yr of preeruption inflation (Figure 6(c)) and 7 cm/yr of coeruption deflation (Figure 6(d)).Figure 6
(a) and (b) are magnified images in Figures4(a) and 4(b) with the cross lines of A, B, C, and D. (c) and (d) represent the profile of averaged surface deformation along each cross line (modified from [16]).


## 4. Volume Change Model of Magma Chamber

The Mogi model, based on a point source embedded in an elastic homogeneous half-space [32], is applied to estimate the location and volume change of the mean surface deformation rate maps (Figures 6(a) and 6(b)) before and during the eruption. The best-fit Mogi source models of the observed interferograms (Figures 7(a) and 7(b)) are represented by synthetic interferograms in Figures 6(c) and 6(d). Residual interferograms (Figures 7(g) and 7(h)) can be calculated from the difference between the observed and synthetic interferograms. The depth of source is estimated to be ~0.3 km below sea level (BSL) before eruption and ~0.6 km BSL during the eruption (Table 2). This Mogi model estimates an increase in volume before the eruption as 1.9 × 10−6 km3/yr (Figure 7). During the deflation of the coeruption period, the volume decrease is estimated to be −2.7 × 10−5 km3/yr (Figure 7). As the Mogi model is essentially a point source, we further generate a best-fit model for using a prolate spheroid [8] for the preeruption and coeruption periods. The prolate spheroid source model (Figure 7(e)) is calculated to be at a depth of 1.3 km BSL beneath the center of the caldera before the eruption and ~1.1 km BSL during the eruption. The RMSE between the observations and model predictions from the prolate spheroid model is less than the corresponding one from the Mogi source (Table 2). Results from the prolate spheroid model suggest that the source of coeruption deflation is a steeply dipping elongate spheroid with major and minor axes of 2.7 km and 1.2 km, respectively, and a loss in pressure of −62.2 MPa/yr (Table 3). The preeruption source axes are estimated to be 3.0 km and 2.4 km with an increased pressure of 4.3 MPa/yr (Table 3).Figure 7
Application of Mogi and prolate spheroid models to the observed mean surface deformation maps during the preeruption and coeruption periods.Table 2
Source parameters obtained from the Mogi model.
Depth (km)
ΔV (km3/yr)
RMSE (cm/yr) between line A and B
Preeruption (inflation pattern)
0.3
1.9 × 10−6
0.18
Coeruption (deflation pattern)
0.6
−2.7 × 10−5
1.25Table 3
Source parameters obtained from the prolate spheroid model.
Major axis (km)
Minor axis (km)
Depth (km)
ΔP (MPa/yr)
RMSE (cm/yr) between line A and B
Preeruption (inflation pattern)
3.0
2.4
1.3
4.3
0.09
Coeruption (deflation pattern)
2.7
1.2
1.1
−62.2
0.58


## 5. LAHARZ Processing

The rationale for the delineation of lahar-inundation hazard zones [33, 34] is derived from the scaling analysis of lahar paths and statistical analysis using the LAHARZ program [35–37]. These analyses generate semiempirical equations that predict the inundated valley cross-sectional area (1) and the planimetric area (2) as a function of lahar volume with proportionality coefficients c and C using volcanic event data of 27 paths of lahars at nine volcanoes.
(1)logA=logC+23logV,(2)logB=logc+23logV,where A is the cross-sectional area, B is the planimetric area, c and C are proportionality coefficients, and V is the volume of the lahar. These calibrated predictive equations involve all the information needed to calculate and plot inundation on topographic maps through an empirical formula. Statistical methods from the LAHARZ program result in experience coefficients, C and c, C=0.05 and c=200 ((3) and (4)).
(3)A=0.05V2/3,(4)B=200V2/3.PF inundation maps are generated in this study using a statistically constrained simulation model to estimate potential hazard areas [25] by adapting methodology from the LAHARZ program [35]. The prediction equations for cross-sectional (5) and planimetric (6) areas of PF hazards for Sinabung volcano are calibrated using data from dozens of events at several volcanoes. The equations are expressed with empirical proportionality coefficients (C=0.05 to 0.1, c=35 to 40) resulting from the adaptation of LAHARZ program statistical analyses [35].
(5)A=0.05to0.1V2/3,(6)B=35to40V2/3.The PF inundation map for the 2010 Sinabung eruption is created in this study using the PF modeling algorithm from the modified LAHARZ program with empirical coefficient values of 0.1 forC and 40 for c. Several values for the volume of the pyroclastic flow in (5) and (6) are tested, including the volume change estimate calculated using the Mogi model with deformation observed in the coeruption interferogram (Table 2).The inundation maps for the PF (Figure8) estimated from the LAHARZ program can be compared to the results of the supervised classification created using the posteruptive Landsat 7 ETM+ image from 30 Jul. 2012. The modified LAHARZ program estimates are calculated using the four different values for the PF volume (red, brown, orange, and yellow colors).Figure 8
The pyroclastic flow inundation area (a) estimated from the revised LAHARZ program (PF simulation) with various volumes (20,000, 27,000, 30,000, and 40,000 m3). The black line (b) is the pyroclastic flow inundation area from the supervised classification of the Landsat image acquired on 30 Jul. 2012 (after the 2010 eruption).
(a)
(b)Among the four different values for volume tested using the modified LAHARZ method, the volume of −2.7 × 10−5 km3 (red color) agrees relatively well with the result from the supervised classification delineation for the inundation limits of the PF deposits from the Landsat image. Region B (volume of PF simulation) in Table 4 is also included within approximately 84% of the PF inundation area delineated using the Landsat image. Region A represents a too small PF simulation area compared with the Landsat image result. Regions C and D, however, have common areas of only 78% and 67%, respectively, between the PF simulations and pyroclastic inundation extents based on the Landsat image (Figure 8).Table 4
Comparison of the common region between the pyroclastic flow (PF) simulation and Landsat supervised classification hazard zone in Figure7.
Pyroclastic flow volume (m3)
Orientation (direction)
Common region between Landsat and PF simulation (pixels)
Pyroclastic flow inundation area from PF simulation (pixels)
Proportion of common region (%)
Region A 20,000
E
Underestimated
Underestimated
Underestimated
SE
SSE
S
Total
—
—
—
Region B 27,000
E
252
284
88.73
SE
110
134
82.09
SSE
90
136
66.18
S
130
138
94.20
Total
582
692
84.10
Region C 30,000
E
256
292
87.67
SE
110
140
78.57
SSE
60
142
42.25
S
136
144
94.44
Total
562
718
78.27
Region E 40,000
E
272
356
76.40
SE
110
170
64.71
SSE
60
174
34.48
S
144
180
80.00
Total
586
880
66.59
Region F 50,000
E
272
402
67.66
SE
110
204
53.92
SSE
60
198
30.30
S
144
204
70.59
Total
586
1008
58.13
Region G 60,000
E
272
474
57.38
SE
110
230
47.83
SSE
60
238
25.21
S
144
226
63.72
Total
586
1168
50.17
Region H 70,000
E
272
500
54.40
SE
110
262
41.98
SSE
60
252
23.81
S
144
258
55.81
Total
586
1272
46.07


## 6. Discussion

InSAR images over Sinabung indicate volcanic inflation before the 2010 eruption and deflation during the eruption. Modeling the observed deformation using Mogi and prolate spheroid sources has allowed us to estimate the magma source location, excessive pressure, the shape of the source, and volume changes during the preeruption and coeruption periods. The model results suggest there is a magma storage underneath Sinabung. The reservoir is located ~1 km BSL (~3 km beneath the peak of Sinabung). The magma reservoir is best characterized with a prolate spheroid with the vertical axis of ~2.7 km and horizontal axis of 1.2 km.The location of the shallow magma reservoir matches the observations of seismicity recorded by the local seismic network. The Indonesia’s Centre for Volcanology and Geological Hazard Mitigation (CVGHM) installed four seismic stations (KWR, SKN, SKM, and MRD; green stars in Figure9(a)) around the crater area of Sinabung after the first eruption on 29 Aug. 2010. These stations recorded earthquakes that occurred on September 6 and 7 (red and blue circles in Figure 9(a)). Figure 9(b) displays hypocenters and the depths of earthquakes that occurred before and after another eruption on 7 Sept. 2010 [38]. A magma reservoir at ~1 km BSL under the peak of Sinabung seems to be consistent with the fact that most earthquakes are located shallower than 1 km BSL (Figure 9). Such a magma reservoir is responsible for the observed inflation before the 2010 eruption and the deflation associated with the 2010 eruption.Figure 9
(a) Hypocentral distribution of volcanic earthquakes at Sinabung volcano from 6 to 7 Sept. 2010. Green stars indicate the locations of seismic survey stations. Red and blue dots indicate hypocenters before and after the eruption of 7 Sept. 2010 (from [38]).
(a)
(b)Pyroclastic inundation hazard maps are generated in this study based on a modified LAHARZ program (PF simulation) using the volume change calculated by the Mogi model from the deformation results observed using stacked interferograms and then are compared with the inundation map from the supervised classification method with Landsat imagery. When we generate these pyroclastic inundation hazard maps, the volume change associated with the deflation of the magma chamber after the eruption was found to be the better-fit value (approximately 86% overlap) based on the comparison with the pyroclastic inundation maps from the Landsat images and the other PF simulations. This study has demonstrated that a volcano hazard inundation map can be generated using the PF model [25] by inputting parameters estimated from models of the magma chamber based on InSAR-derived surface deformation measurements. However, our method has some limitations. First, the eruption volume is calculated from the volume change by modeling the volcanic deflation from InSAR. Due to the vascularity, the eruption volume can be different from the deflation volume [8]. Hence, the volume used for LAHARZ modeling represents the lower bound of the volume of eruptive material. Second, we have ignored other eruptive products such as volcanic ash, which can lower our estimate for the volume used in LAHARZ modeling. Third, we have also neglected other thermodynamic effects on the magma chamber in the derivation of a deflation volcano from InSAR-derived deformation measurements [8].


## 7. Conclusions

In this study, we firstly measure the mean surface deformation using stacked interferograms by ALOS/PALSAR data before and after the 2010 eruption of Sinabung volcano on the Island of Sumatra, Indonesia. The mean surface deformation rate maps are used in the Mogi and spheroid models to estimate the depth, volume changes, and dimensions of the magma source before and during the 2010 eruption. The estimated depth of the magma source is around 1 km BSL; the volume change during the coeruption is calculated. The change in volume estimated from the Mogi model is used to generate a PF inundation hazard map using the modified LAHARZ program. We finally have verified the pyroclastic flow inundation area using the inundation map from the supervised classification method based on Landsat 7 ETM+ imagery. The best-fit result is when the volume for the PF simulation is based on the volume change from modeling the coeruption InSAR deformation map. We conclude that the combination of satellite image analysis and PF simulation using the modified LAHARZ program can accomplish volcano monitoring and pyroclastic flow inundation hazard mapping of active volcanoes.



---

*Source: 8217565-2018-03-19.xml*