Method for estimating the end of the deflation initiated in 2014 at Sinabung volcano, Indonesia, under the assumption that the magma behaves as a Bingham fluid

We herein estimated the end of the ongoing deflation at Sinabung volcano, Indonesia, which began in 2014 under the assumption that the magma is incompressible and behaves as a Bingham fluid. At first, we estimated the temporal volume change of the deformation sources for deflation periods from the end of 2013 until December 2016. The deflation rate is decreasing gradually, and the deflation is expected to cease in the future. Next, we obtained the absolute pressure in the reservoir as a function of time by modeling the magma as an incompressible fluid inside a spherical reservoir with a cylindrical vent and applying the estimated volume change function to the model. We then estimated the yield stress of the Sinabung magma to be 0.1–3 MPa from a previously derived linear relationship between the silica content and yield stress of lava and the measured silica content of Sinabung ashfall from 2010 to 2014, which was reported to be 58–60%. The absolute pressure in the reservoir is expected to fall below the yield stress between August and September 2018, which means that the magma would stop migrating from the reservoir toward the summit under the assumption that it behaves as a Bingham fluid. As a result, the deflation of Sinabung is estimated to end between August and September 2018. In comparison with the 1991–1995 deflation at Unzendake, the total deflation volume of Sinabung is comparable, whereas the duration of the deflation of Sinabung is approximately 1 year longer.


Introduction
Sinabung volcano is an andesitic stratovolcano located approximately 40 km northwest of Lake Toba, Sumatra Island, Indonesia (Fig. 1). Figure 2a, b shows the monthly numbers of deep and shallow volcano-tectonic (VT) earthquakes and the monthly numbers of eruptions and pyroclastic flows, respectively, occurring between 2010 and 2016. Between August and September 2010, phreatic eruptions occurred at Sinabung after at least 400 years of dormancy. The number of deep VT earthquakes increased in July 2013, and phreatic eruptions resumed in September 2013. Eruptive activities gradually became violent, and vulcanian eruptions occurred in November 2013 with an increase in the number of shallow VT earthquakes. On December 22, 2013, a lava dome appeared at the summit. Since the appearance of this lava dome, repeated eruptive activity with pyroclastic flows has occurred.
Continuous Global Positioning System (GPS) observation began at Sinabung in February 2011 (Iguchi et al. 2012). The current distribution of GPS stations is shown in Fig. 1b. Figure 2c shows changes in the baseline lengths until December 2016. Ground extension Open Access *Correspondence: hotta@sus.u-toyama.ac.jp 1 Faculty of Sustainable Design, University of Toyama, Toyama 930-8555, Japan Full list of author information is available at the end of the article initiated in June 2013 prior to the appearance of a lava dome and the detection of subsequent ground contraction. In the previous study, we applied a Mogi model (1958) to the deformation. As for the deflation periods starting in December 22, we divided into two periods (periods 3 and 4), as shown in Fig. 2, based on the rate of change of the baseline length. When rapid contraction was detected in period 3 (December 22, 2013-April 20, 2014), a deflation source was identified beneath the eastern flank of Sinabung at a depth of 8.4 km below sea level (bsl). Its volume change was − 20.51 × 10 6 m 3 . When the contraction rate decreased in period 4 (April 11, 2014-January 16, 2016), a deflation source was identified between Sinabung and Sibayak volcanoes at a depth of 12.2 km bsl. Its volume change was − 88.26 × 10 6 m 3 (Fig. 1b). We also estimated temporal volume change until May 2016 with the location fixed at the position of the deformation source obtained for each period and found that the rate of volume decrease tended to decay gradually (Hotta et al. 2017).
In the present study, we estimated the end of the deflation under the assumption that the magma behaves as a Bingham fluid. A Bingham fluid behaves as an elastic body under pressures below its yield stress and as a viscous fluid at higher pressures. To estimate when the deflation will end, it is necessary to obtain the temporal pressure change P(t) in the reservoir and the yield stress σ yield of magma. Under the condition P(t) ≤ σ yield , the magma does not migrate from the reservoir toward the summit, and the time at which occurs is when the deflation terminates.

Observation and data analysis
We used the data from the continuous GPS stations shown in Fig. 1b. However, stations GRKI and SKNL were nonoperational in early 2014 and early 2016, respectively, and observations at station MRDG and KBYK began in July 2015 and December 2016, respectively. In the present study, we focused on the deflation occurring from the end of 2013 until the end of 2016. Therefore, we excluded the data from stations GRKI and KBYK, and considered only data from stations MRDG and SKNL collected after July 2015 and until January 2016, respectively. We used the L1 and L2 GPS frequency bands and the daily positions with respect to the international terrestrial reference frame (ITRF) 2008 (Altamimi et al. 2011) were calculated by precise point positioning with ambiguity resolution (PPP-AR) analysis (Zumberge et al. 1997) using the software GIPSY-OASIS II ver. 6.1.2. The sampling intervals of GPS data recorded were every 1 s. 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 highprecision ephemerides obtained by the Jet Propulsion Laboratory (JPL) in the present calculations.

Temporal volume change of deflation periods
At first, we estimated the temporal volume change for the deflation from the beginning of period 3 until December 2016 fixing the location of the deformation source obtained by applying Mogi model for each period. We assumed that the location of the deformation source after January 2016 was the same as the source location obtained for period 4. We corrected the topographical effect by adding the altitudes of the stations to the source depth (Williams and Wadge 1998). We divided entire period 3 into two sessions, and the time from the  (Hotta et al. 2017) beginning of period 4 until December 2016 into 20 sessions, respectively. We selected the time window for these sessions such that the most data were available for each, because GPS data were sometimes lost due to transmission problems or power outages as seen in   (Hotta et al. 2017) errors (1σ) were calculated as the square root of the sum of the variance of the first and last 10 days of the time window at each station. We excluded the data from station LKWR obtained during the periods from March to May and August to November 2015 because of a local effect around the station (Fig. 2c). We searched the optimal volume change of the deformation source for each time window which minimizes the residual function f, which is defined as where dx i,j , dy i,j , and dh i,j are the relative east-west, north-south and up-down 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. Figure 3a shows the estimated cumulative volume change of the deformation sources from the beginning of period 3 until December 2016. The observed deformation can be accurately explained by the estimated volume change of the source, as shown in Fig. 3b. The cumulative volume change from the middle of 2015 until December 2016 remained greater than that predicted by the linear trend based on the cumulative volume change from the beginning of period 3 until the end of 2014 (gray dashed line in Fig. 3a), and the rate of deflation tended to decline. In addition, discharge rate of Sinabung is decaying gradually (Nakada et al. 2017) which implies that deflation rate should be also decaying gradually. The deflation is expected to decrease gradually and ultimately cease.

Absolute pressure in the reservoir as a function of time
Next, we obtained absolute pressure in the reservoir P(t) as a function of time using the model by Nishimura (1998) with a spherical magma reservoir , and a cylindrical vent applying Bernoulli's equation.
As Nishimura (1998) assumed, magma behaves as an incompressible fluid at depth. If the magma flow is assumed to be incompressible and the size and location of the reservoir remain constant, the pressure change as a function of time can be written as: where P 0 is the initial pressure in the reservoir, ρ is the density of magma, K is the bulk modulus of magma, S is the cross-sectional area of the vent, and V is the volume of the magma reservoir. Under the condition that the radius of the reservoir is sufficiently small compared to its depth, the magma reservoir can be approximated as a Mogi model. Assuming Poisson's ratio to be 0.25, the volume change as a function of time can be written as (Delaney and McTigue 1994) (3) Fig. 3 a Cumulative volume change of the deformation source from the beginning of period 3 until December 2016 (black dots). The error bars are 95% confidence intervals estimated from the F-test (Árnadóttir and Segall 1994). The gray dashed line shows the linear trend of the cumulative volume change from the beginning of period 3 until the end of 2014. b Comparison of observed (black dots) and calculated (gray dots with dashed lines) changes in baseline lengths where μ is the modulus of rigidity of magma and C is a constant. The total decrease in volume (3VP 0 /4μ) should be sufficiently small compared to the volume V of the magma reservoir because of the assumption in the model by Nishimura (1998) that the size of the reservoir does not change. From Eqs. (3) and (4), temporal volume change can be written as for an andesitic magma (Malfait et al. 2014), and ρ = 2500 [kg/m 3 ], we approximated the cumulative volume change from the beginning of period 4 until December 2016, when the location of the deformation source was assumed to be stable at the source location obtained for period 4, by Eq. (5) using a least-square method. Here, t was set to be the number of days from April 16, 2014. This yielded P 0 = 315 [MPa], S = 707 [m 2 ], and V = 1.58 × 10 10 [m 3 ]. Equation (5) can be fitted to the cumulative volume change, as shown in Fig. 4. Setting the volume of the reservoir 1.58 × 10 10 m 3 yields a reservoir radius of 1.56 km. This radius is less than 20% of the depth of the deformation source (12.2 km bsl), and thus the approximation of the reservoir using a Mogi model is valid (Lisowski 2007). The total volume decrease 3VP 0 /4μ was found to be 1.24 × 10 8 m 3 , which is less than 1% of the volume of the reservoir V. Therefore, the change in the size of the reservoir is negligibly small relative to the total volume. Substituting the obtained parameters into Eq. (3), the temporal pressure change in the reservoir in units of megapascals is obtained as

Yield stress of magma
We then estimated the yield stress of the magma at Sinabung based on its silica content. Anda et al. (2016) reported that the silica content of ashfall at Sinabung from 2010 to 2014 was 58-60%. Furthermore, Hulme and Fielder (1977) have indicated a linear relationship between silica content and yield stress for some types of lava (Fig. 3 in Hulme and Fielder (1977)). The silica content of ashfall at Sinabung reported by Anda et al. (2016) corresponding to a yield stress range of approximately 0.1-3 × 10 −1 MPa in a relationship given by Hulme and Fielder (1977), indication that the yield stress σ yield of magma at Sinabung can be estimated to be 0.1-3 × 10 −1 MPa.

Estimation of the end of deflation
From the function for the temporal absolute pressure P(t) in the reservoir and the yield stress σ yield of magma at Sinabung obtained as described above, we estimated the end of the ongoing deflation. Figure 5 shows a graph of the pressure change function P(t) given by Eq. (6) plotted on a logarithmic scale. The obtained function P(t) falls below σ yield (0.1-3 × 10 −1 MPa) between August and September 2018. Under the assumption that the magma behaves as a Bingham fluid, the magma is estimated to stop migrating from the reservoir toward the summit during this period, when the pressure in the reservoir is expected to fall below the yield stress of the magma. As a result, the deflation of Sinabung is estimated to end between August and September 2018.

Discussion and conclusions
Under the assumption that magma behaves as a Bingham fluid and is incompressible, it was estimated that the deflation of Sinabung will cease between August (3) (gray dashed line). The error bars are as described in Fig. 3 and September 2018, when the pressure inside the reservoir is expected to fall below the yield stress of the magma. The deflation started at the end of 2013 when a lava dome appeared at the summit. Therefore, the duration of deflation was estimated to be 56-57 months (i.e., approximately 4 3/4 years). From Eq. (5) and the obtained parameters, the volume change as a function of time can be written as: During period 3 (December 22, 2013-April 20, 2014), the total deflation was found to be 20.51 × 10 6 m 3 (Hotta et al. 2017). From this and Eq. (7), the total deflation from December 22, 2013, until the termination of the deflation was calculated to be −20.51 × 10 6 + [(−1.24 × 10 8 + C) − V(0)] ≈ − 1.45 × 10 8 [m 3 ]. We fixed the position and depth of the source at the source location obtained for period 4 with a depth of 12.2 km bsl. If it were made 10% shallower/deeper, cumulative volume change would be about 20% less/ more, respectively. If cumulative volume change were 20% less, P(t) would be obtained as Similarly, if cumulative volume change were 20% more, P(t) would be obtained as:   Both of them would fall below σ yield (0.1-3 × 10 −1 MPa) in August-September 2018, similarly to the case for original depth, and 10% depth change do not affect to the estimation.
The confining pressure along the conduit is not included in the equations by Nishimura (1998). Given that the magma remains within the conduit when the magma stops, we need to take the confining pressure into account. Considering this confining pressure, the pressure change in the reservoir as a function of time is re-written as P(t) + ρgh, where ρ, g and h are the density of magma, the gravitational acceleration and the height of magma column, respectively. When ρgh > σ yield , the magma stops without falling below the yield stress of the magma. In this case, the magma stops at t = 1630, i.e., October 2, 2018. Therefore, considering the confining pressure or not do not significantly affect to the estimation of the end of deflation.
The activity of Unzendake volcano, western Kyushu, Japan, from 1991 to 1995 is a case that is similar to the current Sinabung activity; in the Unzendake case, inflation and subsequent gradually decaying deflation occurred after the appearance of a lava dome in May 1991, which was detected from GPS campaign data (Nishi et al. 1999). The deflation of Unzendake is considered to have started in May 1991 when the lava dome appeared at the summit, and the deflation terminated around February 1995. The duration of the deflation was 45 months, which was approximately 1 year shorter than that predicted in the Sinabung case. Nishi et al. (1999) estimated the deflation volume to be approximately 1.5 × 10 8 m 3 , which was almost the same as that for the Sinabung case. This means that the deflation rate during the deflation episodes in the Sinabung case is smaller than that in the Unzendake case.
SO 2 emission has been actually detected at Sinabung ). SO 2 /H 2 O molar ratio at typical volcanoes is about 0.01-0.03 (Shinohara, personal communication on January 12, 2015). These mean that H 2 O is also included in Sinabung magma. If Sinabung magma is with 4-5 weight percent water, the magma would probably gas into bubbles and would be compressible (e.g., Shaw 1974;Neuberg and O'Gorman 2002). Also, we should be aware of CO 2 exsolution which may produce bubble nucleation at even greater depths. However, there is not enough data of CO 2 exsolution to discuss. Hamling et al. (2016) modeled post eruptive subsidence detected from Interferometric Synthetic Aperture Radar during two phreatic eruptions occurred at Tongariro volcano, New Zealand, in 2012. They obtained a shallow deflation source at a depth of 1100 m above sea level. They suggested that fractures associated with the eruptions caused depressurization at the shallow hydrothermal system and subsidence will continue until the fractures become released. Shallow depressurization might also occurred during the eruption and deflation series at Sinabung although such a shallow deflation cannot be detected from our current GPS observation network.
The end of the deflation at Sinabung was estimated from continuous GPS data collected until December 2016 based on the assumption that a new magma intrusion would not occur after the end date of the collected data. On January 16, 2017, earthquakes swarmed near the deformation source obtained for period 4 located between Sinabung and Sibayak, and a change in the ground deformation rate was detected from the automatic analysis of Leica GPS Spider. Ongoing GPS observation and analysis at Sinabung is recommended so that changes in the style of magma intrusion and migration in future activity can be monitored.