Magmatic inflation in 2008 2010 at Mt. Fuji, Japan, inferred from sparsity-promoting L1 inversion of GNSS data

Mt. Fuji, the tallest volcano in Japan, has not erupted since the Plinian eruption in 1707, but its volcanic activity has not ﬁ nished. Here, using GNSS(Global Navigation Satellite System)time-series data, we investigate anin ﬂ ation event of Mt. Fuji in 2008 – 2010. Referring to recent studies of igneous processes, we focus on “ vertically-ex- tensive ” magma systems and estimate depth distribution of geodetic in ﬂ ation sources. A sparsity-promoting L1 regularizationalgorithmforinversionanalysisallowsustoimproveperformanceofdata ﬁ tting.Ourresultshows that the in ﬂ ation sources were mainly distributed in a depth range from 20 km to 15 km, and the total volume change was on the order of 0.01 cubic kilometers (approximately a sphere of radius 100 m). The 15 km depth maycorrespondtotheupperboundaryofpartialmeltofbasalticmagma.Oursparsity-promotingmethodisuse- ful to constrain magma movement from limited observation data.

Nowadays, a dense GNSS (Global Navigation Satellite System) array, named GEONET, by GSI (Geospatial Information Authority of Japan) monitors surface deformation accurately. On the basis of GNSS timeseries data, GSI (2011) tentatively reported an inflation event of Mt. Fuji in 2008-2010 with one point source assuming a fixed depth. Harada et al. (2010) also reported observation of dilatational strain change and discussed its relation to the low-frequency earthquakes.
After that, on March 15, 2011, a large earthquake (M6.4) whose hypocenter depth was about 15 km occurred in the southern foot of Mt. Fuji (Fig. 1). This strike-slip earthquake (hereinafter called as Shizuoka earthquake) occurred 4 days after the 2011 Tohoku megathrust earthquake (M9.0). This Shizuoka earthquake might be related to the static and dynamic stress changes by the Tohoku megathrust earthquake (Fujita et al., 2013;Namiki et al., 2016), and local stressing due to magma activity beneath Mt. Fuji (Hosono et al., 2016).
Recent studies of igneous processes have revealed "vertically-extensive" magma systems beneath volcanoes (e.g. Cashman et al., 2017). Hence depth distribution of geodetic inflation sources must be important. Here, in order to examine the magma activity of Mt. Fuji, we estimate the source of the 2008-2010 inflation of Mt. Fuji by a new analysis method.

Material and methods
We use displacement time-series data of GNSS, the F3 solution of GEONET (Nakagawa et al., 2009), for 15 stations (Table 1) surrounding Mt. Fuji (Fig. 1). We checked precision of three-component (North-South, East-West, and Up-Down) data in the F3 solution, and confirmed that the Up-Down component is approximately three times worse than the other components. Therefore, only horizontal components are analysed in this study. In addition, we did not use a station around the summit of Mt. Fuji (021100), because its time-series data have many missing values.
We first remove offsets owing to antenna replacement and coseismic deformation due to a large event of the Suruga-Bay earthquake (Aoi et al., 2010) on August 11, 2009. Fig. 2a exhibits the timeseries data in horizontal direction away from the summit (138.727°E, 35.361°N) at representative stations. Then the data were detrended for linear, annual, and semiannual components using a prior period as steady state. We selected the prior period from February 2005 (after the end of a slow slip event in an adjacent area (Ochi and Kato, 2013)) to January 2008. Fig. 2b shows the detrended time-series data. The data indicate that the stations tended to move away from the summit approximately in 2008-2010, which implies the inflation event of Mt. Fuji.
In order to identify the starting time of the event, we perform an analysis of piece-wise linear regression (Muggeo, 2003) for estimating a breakpoint of each detrended time-series data. This method allows us to find an inflection point of two linear regression lines by the maximum likelihood estimation. We assume the inflection points during a period from July 2007 to June 2009, and estimate them as a reference of the beginning of the inflation event. Fig. 3 shows the result. On the basis of the regression results, we hypothesize three event periods for comparison, where the length of the period is fixed to be one-and-half years. The most probable model is "Period I" from August 2008 to February 2010 that was presented in Fig. 2b. In a later discussion, we additionally consider half-year later and half-year earlier models, named as "Period I-a" and "Period I-b".
Then we obtain transient deformation amount at each GNSS station, as difference between median values of the first month and the last month. The estimation errors are given by sums of interquartile ranges of the first month and the last month. Fig. 4 shows the transient deformation amounts during "Period I" around Mt. Fuji. The radiation-like pattern evidently indicates an inflation event of Mt. Fuji.
Using the transient deformation amounts, we perform inversion analysis of volume changes in geodetic inflation sources for "vertically-extensive" magma system (Cashman et al., 2017). Since a recent seismological study using receiver function analysis (Kinoshita et al., 2015) did not find sill-like horizontal structures and that using seismic tomography (Nakamichi et al., 2007) did not reveal an entire image of the magma system, we consider that the magma system of Mt. Fuji does not contain well-developed magmatic lenses. Therefore, as a conceptual model for "vertically-extensive" system beneath Mt. Fuji, we set numerical grids of the point sources within an elastic half-space (Mogi, 1958) at depths from 5 km to 30 km with a 5 km interval. The Poisson's ratio is assumed to be 0.25. The horizontal positions of the numerical grids are fixed beneath the summit of Mt. Fuji, because of the limitation of the observation points (we will discuss this point in a later section of 4.1).
We should select an inversion algorithm. Our problem of estimating the inflation sources can be described as follows: the observed transient deformation, d, the volume change of point sources, m, and Green's function within an elastic half-space (Mogi, 1958), G, form a forward relation Gm = d. The well-known least-square solution is m = (G T G) −1-G T d. We first tried this traditional least-square solution. However, the estimated volume changes strongly vary in depth (Table 2), and thus Fuji (red triangle) in Japan. Two oceanic plates subduct beneath the main islands of Japan. The broken lines represent the convergent boundaries (e.g. the Nankai Trough or the Japan Trench). The white squares along the Nankai Trough illustrate the source areas of the 1707 M8.7 Hoei earthquake (Chesley et al., 2012) and those along the Japan Trench illustrate the source areas of the 2011 M9.0 Tohoku earthquake (Nishimura et al., 2011). An enlarged map is shown within the red square. For this part, we use the digital elevation model by Geospatial Information Authority of Japan. The spatial grid size is 10 m. The black points are the GNSS stations of GSI used in this study. The numbers show the station IDs for representative points. The focal sphere indicates the centroid moment tensor of the 2011 M6.4 Shizuoka earthquake, estimated by Japan Meteorological Agency. we judged them as unphysical results. Since d and G include observation noise and modeling error, we need a regularization method to obtain proper solution of m (e.g. Aster et al., 2012). One of the most popular method is the L2 regularization algorithm (damped least-square algorithm) as min ‖Gm − d‖ 2 /2N + λ‖m‖ 2 /2, where N is the observation number (2 horizontal components times 15 stations) and λ is a hyper parameter that controls the weight of the regularization. This algorithm (Friedman et al., 2010) constrains the absolute values of the estimated parameters, which can smooth distribution of the solution. Another method is the L1 regularization algorithm as min ‖Gm − d‖ 2 /2N + λ‖m‖, which can promote sparsity of the solution reducing explanatory parameters (Tibshirani, 1996). This type of the regularization has been introduced in geodetic study for earthquake source estimation (e.g. Evans and Meade, 2012).
In this study, we adopt both algorithms and compare the results in terms of data fitting performance. The hyper parameter λ for each case is chosen by the leave-one-out cross validation (Mosteller and Tukey, 1977), which uses one observation as the validation data, and the other observations as the training data. Table 2 shows the estimated volume changes by the L1 inversion (with the best λ = 0.000195) and L2 inversion (with the best λ = 0.00186), in comparison with the traditional least-square algorithm. We propose that the sparsity-promoting L1 inversion is the best algorithm from the following two reasons. First, the sparsity-promoting L1 inversion better constrains the depth distribution that seems physically reasonable and is consistent with previous studies as discussed in a later   section of 4.2. Second, the result of the sparsity-promoting L1 inversion is statistically better than the others in terms of the Bayesian information criterion (BIC), introduced by Schwarz (1978). Namely, the BIC value for the L1 inversion result is −275.67 (smallest), that for the L2 inversion result is −261.05, and that for the traditional least-square solution is −271.13. This result does not conflict that the smallest meanresidual model is the traditional least-square solution (2.40 mm in horizontal components) in comparison with the L1 inversion result (2.66 mm) and the L2 inversion result (2.74 mm), since BIC penalizes the number of model parameters. Fig. 5 illustrates the calculated result of the Mt. Fuji inflation event based on the L1 inversion (the original data is provided by Supplementary Text S1 and S2). As shown by Table 2, the inflation sources are located in a depth range from 20 km to 15 km, and the total volume change is around 2.1 × 10 7 m 3 (0.021 cubic kilometers). The total volume change of the sources approximately corresponds to ejecta volume of a VEI 3 eruption, which is two orders of magnitude smaller than that of the Hoei eruption (VEI 5), although the volume change of the point source and expected ejection volume would be slightly different due to magma compressibility (Johnson et al., 2000). In addition, our results did not reflect several complex effects such as shapes of the inflation sources (Davis, 1986;Amoruso and Crescentini, 2013), viscoelasticity (Dragoni and Magnanensi, 1989;Bonafede and Ferrari, 2009), or crack opening (Okada, 1985;Fialko et al., 2001). Those effects will be implemented in our future inversion analyses.

Horizontal position of the numerical grids
We fixed the horizontal positions of the "vertically extensive" point sources beneath the summit of Mt. Fuji (138.727°E,35.361°N). In order to validate this assumption, we test other horizontal positions of the "vertically extensive" sources: lattice points for 0.08 and 0.16°from the summit to the north, south, east and west. In each case of the horizontal positions, the L1 inversion is conducted. Fig. 6 shows the results of the mean residuals for all cases (the original data is provided by Supplementary Text S3). We confirm that the original horizontal position (beneath the summit) leads to the minimum mean residual (2.66 mm in horizontal components), and thus it is the best model.

Temporal evolution of the inflation event
Next, we examine the identified inflation period, Period I, from August 2008 to February 2010. We test other possible time windows of a half-year later and earlier period: Period I-a (from February 2009 to August 2010), and Period I-b (from February 2008 to August 2009). For both periods, the sparsity-promoting L1 inversion lead to smaller mean residuals than the L2 inversion. Table 3 summarizes the L1 inversion results of the depth distribution of the volume changes. For the Period I-b, the inflation sources are distributed at depths of 25-20 km, and the total volume changes are about 1.5 × 10 7 m 3 . For Period I-a, the inflation source almost concentrates at a depth of 15 km and the volume change is around 1.2 × 10 7 m 3 . Fig. 7 illustrates the above results.
We can interpret the above results as a simple model that a clot of magma moved upward from a deep part (~25 km, possibly representing the magma reservoir depth) to a depth of 15 km. The depth of 15 km corresponds to the upper boundary of partial melt of basaltic magma (Nakamichi et al., 2007) and the source region of low-frequency earthquakes beneath Mt. Fuji (Nakamichi et al., 2004). A magnetotelluric study also revealed a conductive body at depths greater than 15 km (Aizawa et al., 2004). Furthermore, the hypocenter and lower boundary of the source fault of the 2011 Shizuoka earthquake was also around a Table 2 Inversion results of the depth distribution of the volume changes (the positive values mean inflation) by the traditional least-square solution, with the sparsity-promoting L1 regularization, and ordinary L2 regularization, for Period I from August 2008 to February 2010.
Whether the magma clot stayed or went back after 2010 is an interesting issue. It may be related to the occurrence of the moderate 2011 Shizuoka earthquake. Stress modeling would be useful, however, detailed analyses are rather difficult due to large postseismic deformation by the 2011 Tohoku megathrust earthquake. Development of methods to exclude such high noise is necessary for further studies in this region.

Conclusion
We estimated the depth distribution of the inflation sources beneath Mt. Fuji during 2008-2010, on the basis of the sparsity-promoting L1 inversion of the GNSS time-series data. The estimated sources were mainly distributed in a depth range from 20 km to 15 km, and the total volume change was around 2.1 × 10 7 m 3 . Our sparsitypromoting method is useful to constrain magma movement from limited observation data.