Imaging Volcano Eruption Column Using GPS Tomography: the August 2000 Miyakejima Volcano Eruption
DDudy Darmawan1,2), Fumiaki Kimata2), Kazuro Hirahara2), Bambang Setyadji1) 1) Department of Geodetic Engineering, Institute of Technology Bandung, Indonesia 2) Graduate School of Environmental Studies, Nagoya University, Japan E-mail: dudy@gd.itb.c.id
Abstract
Miyakejima volcano island is an active basaltic andesite stratovolcano located about 200 km south of Tokyo, in the Izu-Bonin volcanic chain. The recent intense seismic activities occurred during the period from June 26 to September 4, 2000, which were accompanied by some eruption events. The largest eruption occurred on August 18, 17:00-19:00 JST, producing an eruption column higher than 10 km. The total amount of volcanic ash and lapilli ejected by such eruption was about 8x109 kg, covering the whole area of the volcano. The aim of this research is to study the eruption collumn of the August 18, 2000 eruption using residual GPS data. GIPSY-OASIS II package was employed to process GPS data. Epoch by epoch analysis of the zero-differenced residual data clearly shows that the signal’s travel was delayed. GPS measurements on August 18 indicate that during 2-hour eruption the recorded GPS signals seemed to have been significantly disturbed by the ejected volcanic materials, for especially the signals that traveled across the caldera. The presence of hot volcanic materials inside the eruption column is mainly considered to be the cause of these delays. According to the residual increment, the delay exceeds as large as 15 cm in some cases. In contrast, the residual decreases (negative) if the signal did not travel across the caldera. The result from residual analysis is then further analyzed to model the eruption column. In this analysis, we apply GPS tomography technique to create the images of the eruption column. The result generally shows that evolution of the intensity of the volcanic material can be reconstructed, which increase during the eruption. According to the reconstructed images, the ejected volcanic materials were likely to be directed to the northwestern of the volcano. Unfortunately due to lack of the inversion model, the eruption column cannot be reconstructed well. 1. Introduction Tomographic techniques aim to obtain a solution field from integrated measurements along different ray paths. They have been successfully applied, for instance, in seismic and oceanographic studies, using earthquake and acoustic sensing; two-dimensional images of the ionosphere using radio signals from polar orbiting satellites and dedicated receivers have also been described. The Global Positioning System (GPS) globally and continuously transmits radio signals, which are affected by the presence of the atmosphere and hence carry information on its state. Thus, GPS signals can be used as data for 4D atmospheric tomography (Flores et al, 2000). Application of GPS tomography in assessing atmospheric studies has been proposed during the past few years. While global ionospheric maps have been described earlier (Ruffini et al., 1998), Seko et al (2000) described regional troposphere four-dimensional maps in Kyushu, Japan, and local scale tropospheric tomography has also been described in several reports (Hirahara, 2000; Flores et al, 2000; Flores et al, 2001; Gradinarsky, 2002). There are many factors that can affect the global atmosphere, one of which is volcanic explosion. Eruption of a volcano may give significant effect to the global atmosphere, as it had been shown by the two-old successive volcanic explosion, Tambora (1815) and Krakatau (1883) (Mills, M.J., 2000). A plinian eruption of Pinatubo (1991) in Philippine also affected the regional atmosphere. In this case, a volcano eruption then provides a unique opportunity for GPS tomography, since it might possible to reveal out 4D maps of the state of eruption column. This research is attempted to study the eruption column of the August 18, 2000 Miyakejima eruption by using GPS tomography. Eruption in August 18 was the largest one among the successive eruptions occurred during the 2000 Miyakejima eruption. 2. Overview of the Miyakejima 2000 Eruption Miyakejima volcano is an active basaltic andesite stratovolcano located about 200 km south of Tokyo, in the Izu-Bonin volcanic chain (Figure 1). The recent intense seismic activity of the volcano occurred during the period from June 26 to September 4, 2000, which was accompanied by a large seismic swarm (June 26 - July 8, 2000), caldera collapse (July 8, 2000) and some eruption events (July 8, 10, 14, 15, August 18, 29) (JMA, 2000). A large ground deformation during the seismic swarm, which corresponds to a dyke intrusion mostly offshore west of the volcano, has been detected (Yamaoka et al, 2000, Irwan, 2003). The Miyakejima eruption in 2000 was characterized by the formation of a collapsed caldera, which was ~1.6 km across and 6x108 m3 in volume in summit area. Intermittent phreatic and phreatomagnetic eruptions from the caldera contrasted with historical volcanic activity. The total mass of eruptive materials was 1.7x1010 kg and mostly consisted of accessory and accidental materials. The eruption was divided into four stages (Nakada et al (2001) as quoted by Geshi et al, 2002): (1) intrusion stage (26 June-8 July), (2) summit subsiding stage (8 July-middle of August), (3) explosive stage (middle of August-September), (4) degassing stage (September). A collapsed caldera grew in the summit area of the volcano between July 8 and the middle of August. The first summit eruption began on July 8 and continued for 10 min with the generation of whitish eruptive clouds about 800 m above the summit. About 1x108 kg of volcanic ash and lapilli covered the eastern part of the volcano. A collapsed caldera, 1 km across, then appeared after the eruption. The second summit eruption occurred in the collapsed crater during the morning of July 14 and continued intermittently till noon of July 15. Dark-colored eruption clouds rose higher than 1 km above the summit. Volcanic ash (~3x109 kg) was erupted and covered the northeastern part of the volcano. The collapsed caldera enlarged in depth and diameter until the middle of August. Phreatic and phreatomagmatic eruptions with ash emission restarted on August 10 and continued until late September to produce about 7x1010 kg of eruptive materials. A relatively large eruption during the afternoon of August 18 produced about 8x109 kg of volcanic ash and lapilli, which then covered the whole area of the volcano. The duration of such an event was relatively long; it was began around 17.00 JST and declined before 19.00 JST. The eruption column rose higher than 10 km (14 km according to weather-radar observation) and drifted southward. Some of the ash cloud reached Hachijojima Island, located 110 km south of the volcano. Ballistic Fragment larger than 1 m across reached the northwestern coast of the volcano. Basaltic volcanic bombs were issued towards the end of the eruption. Eruptive materials contained about 40% fresh basaltic fragments and bombs, which are considered to be juvenile material. According to visual observation at Ako tide gauge station (west coast of the volcano), the sky was cloudy for a few hours after eruption and the maximum visual range was about 3 km. No significant deformation was observed during this event. 3. Data Processing and GPS Tomography Approach 3.1. GPS Data Processing Double frequency GPS data collected from 8 GPS stations around Miyakejima volcano were processed by employing GIPSY-OASIS II software package (Webb and Zumberge, 1993). Those stations are operated by several institutions, 4 stations of GSI (Geographical Survey Institute), 2 stations of JHD (Japan Hydrographic Department) and the rest of ERI (Earthquake Research Institute). GPS data processing is described as follows. As a preliminary step, an accurate position of each GPS station was estimated by processing long term GPS data: the beginning of the explosive stage (August 1 - 17, 2000). The Precise Point Positioning (PPP) module was employed to estimate the position. We used a cut-off elevation angle of 150 and estimation batch interval of 5 minutes. GIPSY uses Kalman filtering technique to model the zenith tropospheric delay (ZTD) and other time-dependent parameters. The ZTD was modeled as random walk with the process noise of 0.17 mm/s1/2, which allows 1 cm2 variation of tropospheric delay after 1 hour. Niell's mapping function (Niell, 1996) was used to estimate ZTD. The tropospheric delay gradient was also estimated. The maximum phase residual allowed in outlier rejection criteria was set to the default value (5 cm), as it is suggested by GIPSY. The resulted position of the stations may be relatively free from the GPS errors. In this case, effects of a small deformation during the beginning of the explosive stage may remain in the position estimates. In the next step, the data on August 18, 2000 were processed by using the same module like in the previous step. Rapid static strategy with 30 seconds estimation interval is used. Only the ZTD was estimated in this case. In order to preserve realistic residual due to unmodeled effects of the volcanic materials, the maximum phase residual allowed in outlier rejection criteria was set to be 20 cm. GPS site positions were held fixed in this processing step. Furthermore, the postfit LC (ionospheric free linear combination) slant phase residual of each GPS satellite data was computed. Using the scenario described above, the residue may reflect the presence of the volcanic materials, and hence it can carry information on the state of the eruption column. These residues were then used as data for the tomographic approach (see Fig. 2). 3.2. Inversion of Tomographic Data In this section, tomographic inversion method applied for an erupted volcano is described. The relation between the computed slant-data residual (Lsr) and the refractivity due to the volcanic materials (Nv) along the slant path is given by (Flores et al., 2001);  If we divide the volume space above the GPS site into 3D finite volume pixels (henceforth voxels), having constant refractivity, eq.(1) transforms to:  where a(i,j) is the length of ray i in voxel j, x(j) is the voxel j refractivity in units mm/km. By exploiting all residual slant observations y for a given time epoch from all available GPS sites and mapping them to the refractivity parameter space x by using the design matrix A, one can form the linear system of equations:  Depending on the used grid the system equations (3) could be ill conditioned due to poor voxel coverage in some areas of the grid defined by the present GPS satellite locations. Some voxels would be over while others under determined. The Square Root Information Filter (SRIF), which is a modified Kalman filter, is applied in the inversion procedure to estimate the refractive index. The SRIF data equations are (Lichten, S.M, 1990):  R, the information matrix, is the inverse of the square root of the covariance. The vector z contains the information needed to obtain the parameter estimates ( ) after inversion of the matrix. ^x ^R We break the data flow into several batches (b) and add every a new batch data (jth) into the previous one (j-1th). We then combine them into an augmented matrix form. By applying the householder transformation to the augmented matrix, we can estimate and for the ^ j R ^ j z  ö is the householder transformation and e is the residual. and values become the a priori information for the next batch. At the end of each batch, Gaussian white noise is added to the parameter estimates ( ), while the random walk process noise (qk) is added to the covariance matrix (P). ^ j R ^ j z ^x 4. Result and Analysis 4.1. Residual of GPS Data An epoch by epoch value of the postfit phase LC residual during the eruption is shown in Figure 2. The residual is plotted 1 hour before, during and after the eruption. Number of the satellites recorded by the stations was 10 satellites. The figure shows clearly that the general residual pattern during the eruption (17:00 -18:00 JST) was noisy and rapidly fluctuated. The fluctuation ceased at 18:30 JST, which corresponds to the end of the eruption at 18:00 JST, followed by some degassing events afterward. Normal value for the residues before the eruption was ± 2 cm, while during the eruption it exceeds 15 cm in some cases. It is of interest that the residue has a large positive and negative value during the eruption. Furthermore, we analyze the cause of such a residual pattern by examining the satellite's travel when it was being recorded by the station. According to the report issued by VRC (Volcano research Center)-ERI (Earthquake Research Institute), the abundant ash felt mainly in the northwestern part of the volcano island, indicating that the volcanic materials were ejected toward the northwest. We then analyze a satellite (SV5) observed at the eastern part of the caldera and two satellites (SV21 & SV 23) at the western part (Figure 3a). SV5 orbited from the south to the northeast while other satellites orbited from the north to the south. Besides their relative positions to the caldera, those satellites are chosen due also to their high elevation angles, above 30o, during the eruption day. Figures 3b, 3c and 3d show the residual time series of SV5, SV23 and SV 21 respectively, from 17:00 JST until 19:00 JST, recorded at station 3059, 0599, 3060 and 0600. The residual value for SV5 in 0600 (Figure 3b) has the largest one compared to other sites. The residual for this site looks very noisy, where the residual becomes large around 17:55 JST. Maximum fluctuation of the residuals is about 15 cm. The same large value also appears in the residual for SV23 in 3059 (Figure 3c) and for SV21 in 0599 and also in 3060 (Figure 3d). Although pattern of the noisy residuals is little complicated, but it shows general features. Gradually the residuals increase after 17:12 JST, except for SV23 in 3059 increase after 17:35 JST, and attain their maximum values around 17:50-18:00 JST. The time interval of the increment of the residuals is consistent with the time interval of the eruption. It may indicate that the erupted materials have given a significant effect on GPS data. To make it clearer, we examine the directions of the received GPS signals as shown in Figure 3a. The satellites with noisy residuals pass through the eruption column before the signals arrive in the GPS sites. While passing through such a column, the signals are disturbed accordingly. This will be clear if we compare for instance the residual for SV5 in 0600 with 3059. SV5 is observed in the eastern part of Miyakejima (azimuth 45o to 165o) and its signal travels to 3059 without passing through the eruption column, but, in contrast, it should pass through the eruption column before arriving in 0600. Therefore, as shown in Figure 3b, the residual in 3059 is much less noisier than 0600. Magnitude of the fluctuation is about 5 cm, three times lower than for 0600. On the other side, the residual will be negative (decrease) if the signal does not pass through the eruption column, as it is shown in Figure 3d for station 3059. The cause of this is being investigated. 4.2. Tomographic Inversion Using the residual generated in the previous section, tomographic inversion is then calculated to determine the refractive index for each voxel above the volcano. Since the height distribution of the station is not good enough to obtain good vertical resolution, the simulations test have been carried out beforehand. In order to determine the impact of the discretisation error, we first define a regular grid of 20 x 20 voxels in the horizontal and 20 layers of 500 m thick up to 10 km height above the lowest GPS antenna to generate the simulated delays. Additive Gaussian noise has been added to the delays. The time period from 16:30-18:00 JST is analyzed by using Kalman filter every 3 minutes (30 batches). For each run, the rms value of the residues has been computed. The result shows that when a coarser grid in the horizontal is used (10 x 10 x 8) the discretisation is below the additive noise. Furthermore, we compute the real data by using the grid generated from the simulations test. In order to account for the time variation of the parameters, a random walk stochastic process with a drift rate of 0.26 cm/sqrt(sec) is used. Obviously, the choice of the rate is a crucial problem in Kalman filter. If the process noise is chosen to small, the Kalman filter will work similar to a low-pass filter. Fine structures in the refractive index will not be visible in the results. Otherwise, if the process noise is chosen too high, the observation will get much weight during the update process. In this case, the Kalman filter does not act any longer as a filter. We choice the process noise rate based on the fluctuation rate of the residual. We assume the process noise rate is the same for the entire Kalman filter batch. Evolution of the refractive index value is shown in figure 4. We plot the value from several batches (8 batches) estimated in the period of 17:06 JST until 17:30 JST. Although the resolution of the inverse process is not so good, we can still recognize that the refractivity of every batch becomes large as the time of batch closes to the eruption time. Several minutes before the eruption (Figure 4a), the refractivity is normal in many places, except for some voxels at the higher layer. The refractivity increases slowly 6 minutes after the eruption (Figure 4b and 4c) and it rapidly increase from 17:15 JST until 17:30 JST. The eruption column cannot be recognized unfortunately. As we can see, when the fluctuation of the intensity is being large, 6 minutes after the eruption, high refractivity mostly occurs in the western and eastern part of the image, while low refractivity occurs near the caldera. According to visual observation of the volcanic plume, the eruption column took place above the caldera. This cannot be reconstructed in the images. According to analyses above, it seems that Kalman filter can not track the time variation of the refractivity. This may due to the process noise rate suppresses the update process for batch interval 6 minutes after the eruption. Therefore, some methods (Segal and Mathews, 1997) to select the optimum process noise rate for each batch should be used in the computation. Although the eruption column cannot be reconstructed, from the period of 17:09 JST until 17:12 JST (Figure 4b), we can still see large refractivity in the northwest of the horizontal image. It may indicate the location of the eruption column occurred at the northwestern part of the volcano as it has been reported by VRC-ERI (http://hakone.eri.u-tokyo.ac.jp/vrc/ erup/miyake.html). 5. Conclusion Using the GIPSY software we processed GPS data during the August, 18 200 Miyakejima eruption. An epoch-by-epoch analysis of the postfit phase LC residual shows that the presence of the hot volcanic material affects GPS signal significantly. The signal pass over the volcanic material seems to have been delayed, since its residual increase rapidly. The delays exceed 15 cm in some cases. In contrast, the residual will be negative (decrease) if the signal does not pass through the eruption column when it is traveling to the stations. The cause of this is being investigated. We then used GPS tomography technique to reconstruct the eruption column. Generally, evolution of the intensity of the volcanic material can be reconstructed, which increase during the eruption. Although the eruption column cannot be reconstructed in the images, we can still recognize the possible eruption column, which is located in the northwestern part of the volcano. Moreover, further studies on inversion procedures needs to be carried out in order to find a reliable inversion model for our GPS data. Our simple inversion model seems not too satisfy as it fails to reconstruct the eruption column. References
    Figure 2. Postfit LC phase residual for all the observed satellites (10 satellites) recorded at all the stations in Miyakejima. Red vertical bar indicates the start of the eruption, which was 17:02 JST according to (JMA, 2000). Red line is the residual for SV5, while cyan for SV9, green for SV30, yellow for SV23, blue for SV29, black for SV21, magenta for SV26 and magenta plus line for other satellites (SV6,17,25).     Figure 3. Sky plot of the satellites and their postfit phase residuals. a) Sky plot for SV5, SV23, and SV21 above Miyakejima volcano from 17.00 JST until 19.00 JST. Solid diamond represents the satellite track in the period of 17:15-18:15 JST. Open rectangular is the location of the abundant ash, reported by VRC-ERI. b) Postfit phase LC residuals for SV5 and c) SV23 d) SV21.     Figure 4. Evolution of the refractive index during the eruption estimated by Kalman filter. a) 17:06-17:09 JST, b) 17:09-17:12 JST, c)17:12-17:15 JST, d) 17:15-17:18 JST, e) 17:18-17:21 JST, f) 17:21-17:24 JST, g) 17:24-17:27 JST, h) 17:27-17:30 JST. The color bar indicates the refractive index (cm/m). Unit for the space model is in kilometer.
Page 1 of 1
|