Title Magma intrusion and effusion at Sinabung volcano, Indonesia, from 2013 to 2016, as revealed by continuous GPS observation

from 2013 to 2016, as revealed by continuous GPS observation Author(s) Hotta, Kohei; Iguchi, Masato; Ohkura, Takahiro; Hendrasto, Muhamad; Gunawan, Hendra; Rosadi, Umar; Kriswati, Estu Citation Journal of Volcanology and Geothermal Research (2017) Issue Date 2017-12-22 URL http://hdl.handle.net/2433/241767 Right © 2018 Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Type Journal Article Textversion publisher


Introduction
Sinabung volcano is an andesitic strato-volcano 2460 m in height. This volcano is located approximately 40 km northwest of Lake Toba, formed by catastrophic caldera eruptions (ca. 74 ka; Chesner and Rose, 1991), in northern Sumatra of Indonesia ( Fig. 1(a)). Lava flows, pyroclastic deposits, and summit domes comprise the edifice of Sinabung volcano (Prambada et al., 2010). The youngest block-andash flow and associated surge deposits, which are distributed on the southeastern flank of the volcano and reach approximately 5 km from the summit, occurred in the 9th and 10th centuries (Iguchi et al., 2012).
In August and September 2010, phreatic eruptions occurred at Sinabung after at least 400 years of dormancy. Ground inflation before the 2010 eruptions and subsequent rapid deflation after the eruptions were detected by interferometric synthetic aperture radar (InSAR) (Chaussard and Amelung, 2012;Chaussard et al., 2013). After the 2010 eruptions, phreatic eruptions started again on September 15, 2013. Eruptive activity increased, becoming phreatomagmatic by November 11, 2013 (Nakada et al., this issue), and the largest vulcanian eruption occurred on November 23, 2013, at which time the alert level was upgraded from three to four (the highest level). On December 18, 2013, lava appeared at the summit of Sinabung (Pallister et al., this issue). Since January 2014, pyroclastic flows have occurred repeatedly as part of the eruptive processes. For detailed descriptions of the eruptive chronology see Gunawan et al. (this issue) and Pallister et al. (this issue). Lava dome/flow and collapse is similar to effusive phases at Sinabung in the 9th and 10th centuries and to the 1991-1995 effusive eruptions at Unzen volcano, located in Shimabara peninsula of western Kyushu, Japan (Nakada et al., 1999).
Continuous Global Positioning System (GPS) observation began at Sinabung in February 2011 (Iguchi et al., 2012). GPS is a useful ground deformation tool to understand magma migration because of its good time resolution, such that it has been successfully used at volcanoes in a number of previous studies (e.g., Sato and Hamaguchi, 2006;Hotta et al., 2016). We use GPS observations at Sinabung to model magma migration and effusion before and during the onset of eruptive activity in September 2013.
In the present study, we analyzed GPS data from stations installed at Sinabung in order to obtain the location and volume change of deformation sources for each deformation stage of recent Sinabung activity since 2013. We then interpreted the obtained deformation sources based on volcanic activity in order to better constrain the magma migration process. Furthermore, we estimated the temporal volume change of the deformation source in order to provide an estimate of the quantitative change in the magma migration and of the termination of effusion. Finally, we compared the current Sinabung activity since 2013 with that during the August-September 2010 eruptions and with that at Unzen volcano during time periods with high rates of pyroclastic flows in 1991-1995.

Observation and data
The GPS network at Sinabung is shown in Fig. 1(b), and the location and data availability for each station are shown in Table 1. While there are five GPS stations at Sinabung, only four or fewer stations have functioned at any given time; station GRKI was nonfunctional in early 2014, and observation at station MRDG only started in July 2015. We used the L1 and L2 GPS frequency bands and calculated daily positions in reference to ITRF2008 (Altamimi et al., 2011) by precise point positioning with ambiguity resolution (PPP-AR) analysis (Zumberge et al., 1997) using GIPSY-OASIS II ver. 6.1.2. The sampling rate of the GPS data is 1 Hz. The GIPSY solutions were run on decimated 5-min intervals, and daily smoothed averages were considered. We used the Global Mapping Function (GMF; Boehm et al., 2006) based on numerical weather model data and high-precision ephemerides obtained by the Jet Propulsion Laboratory (JPL) in the present calculation. Fig. 2(a) shows examples of the temporal change in the baseline length. The ground deformation can be divided into four periods based on observed changes in the rate of ground deformation. In June 2013, slight extension started especially between stations SNBG and LKWR (period 1). The extension accelerated in late October 2013, especially between stations SKNL and LKWR and between stations SNBG and LKWR, with no b1 cm extension in the baseline length (period 2). In January 2014, the baseline lengths between stations SNBG and LKWR and between stations SNBG and SKNL began to contract rapidly at a constant rate until April 2014, resulting in a final contraction of approximately 3 cm (period 3). After April 2014, the rate of contraction of the baseline length between stations SNBG and LKWR and between stations SNBG and SKNL decreased, and the contraction continued until the middle of 2016 (period 4). Changes in the baseline length, which appeared in baselines including the station LKWR, especially during the periods from March to May and August to November 2015, are considered to be a local effect around station LKWR, because these deformations appear only at this station and ultimately returned to the original trend. Another possibility is that the antenna at LKWR may have become coated with ash during these periods, unfortunately signal-to-noise (SNR) was not recorded in our RINEX data therefore we could not test this hypothesis. In addition, during this time period the baseline length between stations SKNL and LKWR, which are located near the summit of Sinabung, has shown only small changes since 2014. Relative tectonic  deformation within our observation network is considered to be negligibly small, because the distance between each station is within approximately 12 km and no clear relative displacement was detected before slight extension began in June 2013. Fig. 2(b) and (c) show the monthly numbers of deep and shallow volcano-tectonic (VT) earthquakes and those of eruptions and pyroclastic flows, respectively, between September 2010 and June 7, 2016. Deep VT earthquakes are high frequency earthquakes ≥ 2 km below sea level (bsl) deep with clear onsets and typical S-P intervals of 1-3.5 s; and shallow VT earthquakes are high frequency earthquakes b2 km bsl with S-P intervals of b 1 s and do not include low-frequency (LF) or hybrid events (Gunawan et al., this issue). After the August-

Methods
Currently, a maximum of four GPS stations are available at any given time at Sinabung. Therefore, it is difficult to consider complex models, such as a prolate spheroid source model (Yang et al., 1988), a rectangular tensile crack model (Okada, 1992), a penny-shaped source model (Fialco et al., 2001), or realistic underground structures (e.g., Hautmann et al., 2010). Therefore we selected to ignore complexities and apply a Mogi model (Mogi, 1958), which is a simple point source in an elastic and homogeneous half-space. This model has been widely used to understand relationships between volcanic activity and magma migration from ground deformation observations (e.g., Nishi et al., 1999;Nakao et al., 2013;Hotta et al., 2016). Of course, we should be aware of the limitations of a Mogi model: for example, the simplified geometry of the deformation source or that it neglects underground structures. We corrected for topographic effects by adding the altitudes of the stations to the source depth (Williams and Wadge, 1998). The ground deformation was divided into four periods, as shown in Fig. 2 and Table 2, based on the rates of ground deformation described in Section 3. We chose the four time windows such that the most data was available for each, because data from one or all stations was lost because of transmission problems and power outages as seen in Fig. 2(a). The observed displacements of each period were calculated as the residual of the average position of the last and first 10 days of the time window at each station; and their errors (1σ) are calculated as the square root of the sum of the variance of first and last 10 days of the time window at each station. The observed displacements and errors for each period in reference to ITRF2008 are as shown in Table 2. We modeled these observed displacements for each period. In order to minimize the difference between observed and calculated values, we defined the residual function f as where dx i,j , dy i,j and dh i,j are the relative east-west, north-south and updown displacements at station i with respect to station j, respectively; n is the number of stations; the subscripts obs and cal represent observed and calculated values, respectively; and σ represents the standard deviation of the observation. σdx i, j , σdy i, j and σdh i, j are calculated as where σdx i , σdx j , σdy i , σdy j , σdh i and σdh j are the standard deviations of the displacements in accordance with ITRF2008 for each component and station. Using a grid search method, we searched optimal parameters of position, depth and volume change of a Mogi source that minimize the function f. The search range and step for each parameter were set as shown in Table 3.

Results
The obtained optimal parameters and 95% confidence intervals for each period are shown in Table 4, and location maps of the obtained deformation source comparing the observations and calculations are shown in Fig. 3. Although there were some differences between the observations and calculations, all of the differences were within  Table 4 Optimal parameters with uncertainties of deformation sources for periods 2-4. The origin of the horizontal coordinate is as described in Table 3. The uncertainties are the 95% confidence intervals estimated from the F-test (Árnadóttir and Segall, 1994 observation errors of 1σ. In addition, the obtained deformation sources are summarized in Fig. 4. We indicated the uncertainties of parameters as the 95% confidence intervals estimated from the F-test (Árnadóttir and Segall, 1994). The uncertainties of horizontal locations and depths are not symmetric about optimal values because their f values do not change symmetrically. For period 1 (June 1 to October 30, 2013), an inflation source was obtained at a depth of 3-8 km bsl with a volume change between + 0.3 and + 1.8 Mm 3 . Although the volume increase at this depth was modeled, the consistency of the position of the obtained inflation source was not good because the observed deformation during this period was small.
For period 2 (October 21 to December 31, 2013), an inflation source was obtained at a depth of 0.9 (0.4-2.1) km bsl directly beneath Sinabung. Its volume change was +0.39 (+0.18-+0.60) Mm 3 . The inflation source in period 2 was shallower than that in period 1.
For period 3 (December 22, 2013 to April 20, 2014), a deflation source was obtained at a depth of 8.4 (7.4-9.9) km bsl beneath the eastern flank of Sinabung. Its volume change was − 20.51 (− 26.89 to − 14.12) Mm 3 . From period 2 to period 3, the volume change transformed from inflation to deflation in this period, and the source shifted eastward and to a greater depth.
For period 4 (April 11, 2014 to January 16, 2016), a deflation source was obtained at a depth of 12.2 (10.1-14.8) km bsl between Sinabung Fig. 3. (a) Obtained deformation source (black dot) for period 2 (October 21 to December 31, 2013). The error bar represents the 95% confidence interval estimated from the F-test (Árnadóttir and Segall, 1994). Arrows indicate the observed (black) and calculated (white)  and Sibayak volcanoes. Its volume change was −88.26 (−123.87 to − 52.66) Mm 3 . From period 3 to period 4, the deflation source shifted northeastward and to a greater depth.

Interpretation of the deformation sources and the construction of the magma migration process
We obtained inflation and deflation sources for the periods of extension and contraction of the baseline length, respectively. The inflation source for period 1 was deeper than that for period 2, which was located at the summit of Sinabung; and the deflation sources for periods 3 and 4 were located east of the summit. Here, each deformation source is interpreted considering relationships between ongoing volcanic processes and the magma migration process as inferred from GPS modeling.
For period 1, an inflation source was obtained at a depth of 3-8 km bsl, but the consistency of the position of the inflation source was not good because the observed deformation during this period was only a slight extension (b 0.5 cm for the change in the baseline lengths). Nevertheless, the estimated error range of the source depth reached a greater depth than that obtained for period 2 (0.4-2.1 km bsl). In addition, the monthly number of shallow VT earthquakes with hypocenter depths of b2 km bsl (Gunawan et al., this issue) was low while deep VT earthquakes were dominant during period 1, in contrast to the subsequent rapid inflation period when shallow VT earthquakes frequently occurred. The inflation of period 1 is therefore considered to correspond to magma intrusion at a location deep beneath Sinabung, precursory to the resumption of phreatic eruptions on September 15, 2013.
During period 2, an inflation source was obtained at a shallow depth (0.9 km bsl) beneath Sinabung. Eruptive activity at Sinabung transistioned from phreatic to phreatomagmatic during this period and included the largest vulcanian eruption on November 23, 2013, and the appearance of a lava dome at the summit on December 18, 2013 (Pallister et al., this issue). In addition, shallow VT earthquakes occurred frequently, especially in November 2013. Furthermore, LF and hybrid events increased dramatically in December 2013 before the appearance of the lava dome (McCausland et al., this issue). This inflation source is therefore considered to correspond to magma migration to shallower depths and the summit from the deeper location where magma was intruded during period 1.
During period 3, the lava dome that appeared on December 18, 2013, grew continually through the movement of magma to the summit. Pyroclastic flows occurred frequently as a result of the collapsing of the lava dome. The obtained deflation source at a depth of 8.4 km bsl beneath the eastern flank of Sinabung is considered to correspond to a magma reservoir, and magma migrated toward the summit from this reservoir during the eruptive episodes occurring in this period.
For period 4, the ground deflation rate was less than that of the previous period. However, pyroclastic flows still occurred frequently. The obtained deflation source was located further northeast and deeper (12.2 km bsl) than that for the previous deflation period. Its depth and volume change magnitude were the greatest among the four deformation stages. The obtained deflation source is considered to correspond to a deeper magma reservoir, which may be the main reservoir for the magmatic supply system of Sinabung.
We tested a finite spherical source (McTigue, 1987) and a vertical spheroid source (Yang et al., 1988;Newman et al., 2006), during period 2, when 4 GPS stations are available and the source depth is very shallow (0.9 km bsl). For a finite spherical source model, the optimal radius obtained is b10 m (horizontal location, depth and volume change are the same as that for Mogi model). In this case, where the radius b 20% of depth, we cannot distinguish a finite spherical source from a Mogi model (Lisowski, 2007). For a vertical spheroid source, the optimal parameters are obtained as below assuming a modulus of rigidity of 30 GPa and a Poisson's ratio of 0.25: horizontal location: E-W 1.1 km, N-S 0.7 km (98.40184°E, 3.17633°N); depth: −0.2 km bsl (i.e., 0.2 km above sea level (asl)); equatorial and polar radii: 0.6 km and 2.3 km, respectively; pressure change: +0.59 MPa. However, the resolved source is unrealisticthe uppermost part of the source is above the ground surface with an altitude of as much as 1.8 km asl at (98.40184°E, 3.17633°N), although one might interpret it as a magma conduit. For these reasons, it is difficult to consider complex source geometries for period 2.
For periods 3 and 4 when only 3 GPS data are available, it is difficult to constraint parameters of complex geometriesfor example, the optimal value of depth and radius of finite spherical source (period 3: 11.0 km bsl depth with 19.1 km radius is unrealistic; period 4: 12.2 km bsl depth with b0.1 km radius cannot be distinguished from a Mogi model)because of the limited station distribution. Fig. 5 shows a schematic of our modeled magma migration process at Sinabung. From June 2013, when slight ground inflation started, magma is considered to have been injected in a location deep beneath Sinabung, followed by an increase in the incidence of deep VT earthquakes ( Fig. 5(a)). The injected magma then migrated to a shallower  Fig. 3, and solid triangles are as described in Fig. 1. location nearer to the summit, and the number of shallow VT earthquakes increased, which caused rapid ground extension at the stations near the summit and finally resulted in the appearance of a lava dome on December 18, 2013 ( Fig. 5(b)). After the appearance of the lava dome, the magma reservoir beneath the eastern flank of Sinabung at a depth of 8.4 km bsl deflated due to magma extrusion and accompanying pyroclastic flows, which caused rapid ground deflation (Fig. 5(c)). From the middle of April 2014, when the ground deflation rate decreased, deflation caused by magma extrusion shifted to a deeper reservoir located at a depth of 12.2 km bsl between Sinabung and Sibayak volcanoes (Fig.  5(d)). This deep magma reservoir may be the main magmatic reservoir for Sinabung.

Temporal volume change of the deformation source
We estimated temporal volume change of the deformation source fixing the location at the modeled position for each of the four periods. We assumed that the location of the deformation source after January 2016 was the same as the source location obtained for period 4. We divided entire periods 2 and 3 into two sections each, and the time from the beginning of period 4 until May 2016 into 16 sections, respectively. The time window for these sections were selected to yield section intervals that were roughly the same for each period (approximately 30 and 50 days on average for period 2 and the time from period 3 until May 2016, respectively); to maximize the available data; and to avoid the use of data at LKWR from March to May and from August to November 2015, when a local effect clearly affected the data at that station. Fig. 6 shows the estimated cumulative volume change of the deformation sources from the beginning of period 2 until May 2016. The observed deformation can be accurately explained by the estimated volume change of the source, as shown in Fig. 7.
During period 2, there was a slight increase in volume until November 2013, and the volume then a greater increase after the vulcanian eruption on November 23, 2013. These results are consist with the magma supply rate toward the summit which was lower during the first half of period 2 and significantly increased after the November 23 vulcanian eruption as the lava made its way to the surface.
During the deflationary sections from period 3, the cumulative volume change from the beginning of period 3 until the end of 2014 followed a linear trend, as shown in Fig. 6. The cumulative volume change after the middle of 2015 was greater than that predicted by the linear trend, and the rate of volume change tended to decrease starting in the middle of 2015. There is some observational evidence that the effusion rate and/or volume change can decrease exponentially (e.g., Wadge, 1981;Lu and Dzurisin, 2010;Anderson and Segall, 2011). Koyaguchi (2016) obtained the following time-dependent functional relationship for the pressure in the reservoir, P, using a model for an over-pressurized spherical magma chamber with a cylindrical magma conduit and assuming the magma flow to be a Poiseuille flow: where P lith is the lithostatic pressure; P 0 is the initial pressure of the reservoir; and τ is a constant. When the magma reservoir is approximated as a Mogi model, the volume change ΔV can be written as a function of the pressure change ΔP as (Delaney and McTigue, 1994) where r is the radius of the magma chamber and μ is the modulus of rigidity. From Eqs. (3) and (4), the temporal volume change can be approximated by a function exponentially decaying with time, as: where V is the cumulative volume change in Mm 3 and t is time in days (for this paper, the number of days after January 1, 2014). We obtained the coefficients in Eq. (5)  ]. Fig. 8 shows a comparison of the estimated cumulative volume change and the estimated exponential function. The exponential function accurately represents the cumulative volume change. Approximately 2/3 of the total volume change related to con-  with a modulus of rigidity of μ = 30 GPa (personal communication). According to Battaglia et al. (2013), the volume change of the spherical source ΔV can be obtained as From Eq. (6), the volume change from 2006 to 2009 was obtained as +0.63 Mm 3 , which corresponds to an average of +0.018 Mm 3 /month. The location of the 2010 inflation source was close to that for period 2 (from October 21 to December 31, 2013), and the average rate of the volume change was approximately 10 times larger in period 2 than that modeled prior to the 2010 eruptions. The characteristics of the inflation and subsequent rapid deflation were similar to those detected by GPS observation during the eruptive activity after September 15, 2013. We do not know if a similar magma migration process occurred during the 2010 eruptions because there was no GPS data available at that time, and the time resolution of the available ground deformation data (InSAR) was poor; importantly no magma, neither juvenile ash nor lava dome, appeared at the surface during the 2010 eruptions. [1991][1992][1993][1994][1995]. Phreatic eruptions at Unzen volcano started in November 1990 and were followed by lava dome formation and subsequent frequent pyroclastic flows beginning in May 1991 (Nakada et al., 1999). This eruption style is similar to the activity of Sinabung starting in September 2013. Nishi et al. (1999) detected inflation of a source from January to March 1991 prior to the dome formation, and obtained a deflation source at a depth of 11 km bsl beneath the western flank of Unzen from GPS campaign data during the lava dome extrusion period from 1991 to 1994. This is very similar to the magma migration and effusion processes of the current Sinabung eruption. Applying Eq. (5) for the cumulative volume change for Unzen estimated by Nishi et al. (1999), we obtained b value as 1.664 × 10 −3 [day −1 ], which is larger than that of Sinabung (b = 1.235 × 10 − 3 [day − 1 ]), which means that the exponential function for Unzen approaches its asymptote more rapidly than Sinabung's. In addition, the maximum decrease in volume estimated from the exponential function for the Unzen case is 150 Mm 3 , which is slightly smaller than that for the Sinabung case (165 Mm 3 ). Therefore, the duration of effusion for the Sinabung activity ought to be longer than that for Unzen 1991Unzen -1994 Assuming the rate of inflation was constant until lava appeared and that the total volume of supplied magma was equal to that of lava discharged, Nishi et al. (1999) estimated that magma intrusion began in December 1989, 17 months before the appearance of the lava dome in May 1991. However, seismic evidence from distal VTs indicates that the intrusions likely began as early as 1984 and occurred in pulses over the 6 years before the onset of phreatic eruptions (White and McCausland, 2016;Nakada et al., 1999;Umakoshi et al., 2001).

Unzen volcano activity in
In the Sinabung case, the increase in the incidence of deep and/or distal VT earthquakes started as early as December 2010 (Fig. 2(b)), with a significant increase in July 2013 (McCausland et al., this issue), 1 month after slight ground inflation started and deep magma intrusion is considered to have begun. In the Sinabung case, the lava discharge volume (~100 Mm 3 ; Nakada et al., this issue) had already exceeded the magma supply volume during the inflation period from June to December 2013 (~2 Mm 3 ). As described in Section 5.3.1.1, magma intrusion should have also occurred before the August-September 2010 eruptions, and the magma intrusion is considered to have already started before 2007 because InSAR data show ground inflation from 2007 at the latest until the start of the 2010 eruptions. The inflation before the 2010 eruptions at Sinabung likely started before 2007 and is a clear sign of magma intrusion before the 2010 and 2013 eruptions. Given the location of the period 4 source and current observation network, inflation rates in the LOS direction detected by InSAR and change rates of baseline length between any two GPS stations are at most 1 cm/year and 0.3 cm/year, respectively, and are difficult to detect within their accuracies when magma intrusion rate is below 7 Mm 3 /year. Fig. 9 shows a schematic of the cumulative volume change resulting from Sinabung activity over time. Assuming that the total amount of magma injection is equal to that of magma ejection, magma intrusion must have started before 2007 as described in Section 5.3.1.2. Between when the GPS network was installed in 2011 and May 2013, there was no measurable deformation at Sinabung, neither regional tectonic nor volcanic. In June 2013, magma intrusion started again leading to the resumption of phreatic eruptions on September 15, 2013, becoming phreatomagmatic in early November 2013 (Nakada et al., this issue), culminating with the largest eruption on November 23, 2013, and finally appearing at the surface as a lava dome on December 18, 2013. Since 2014, magma effusion continued sometimes accompanied by pyroclastic flows. We estimate based on current conditions that magma effusion will terminate in the early 2020s.

Re-consideration of current Sinabung activity and observation
Currently, our GPS observation network is within approximately 12 km around Sinabung volcano, which prevents us from estimating absolute tectonic effects with this network. Although relative tectonic deformations are considered to be negligibly small in our current narrow observation network as mentioned in Section 3, we recommend that the current GPS network be extended in order to estimate absolute tectonic effects around Sinabung so that we may obtain the pure absolute volcanic deformation, which would improve the depths of our modeled sources. The migration of the source to northeast between periods 3 and 4 could be an artifact of not accounting for the overall plate motion. In addition, the current GPS network around Sinabung is still coarse, and thus it is difficult to consider complex models. The Mogi model we used in the present study is a simple point source model and thus unrealistic. We would expect that future deflation would locate similarly to our modeled period 4 source between Sinabung and Sibayak, rather than around Sinabung volcano. Therefore, we recommend that stations be added to constrain the period 4 source as well as the period 1 and 2 sources in order to best model complex source geometries, depending on if the system undergoes continued deflation or has renewed intrusions.

Conclusions
Inflation and deflation sources were detected during periods of baseline length extension and contraction from June to December 2013 and since 2014, respectively. The modeled 2013 inflation periods were similar to the modeled inflation prior to the 2010 eruptions. In 2013, the inflation sources migrated to a shallower location beneath the summit of Sinabung before a lava dome appeared in December 2013. The deflation sources moved eastward and deeper, and began days after the dome appeared at the surface and continue until this writing.
We approximated the temporal change in the cumulative volume during the deflationary periods with a function decaying exponentially with time. The exponential decay function allows the estimation of the total decrease in volume and the estimation of the end of the eruption, assuming no changes in rate or no new intrusions.
The current Sinabung activity is similar to that of Unzen in 1991-1995 in terms of the ground inflation and subsequent gradually decaying deflation with effusive eruption.