Joint inversion of strain and tilt data using the Akaike’s Bayesian information criterion to map detailed slip distributions of short-term slow slip events

Along the Nankai Trough subduction zone, southwest Japan, short-term slow slip events (SSEs) are commonly detected in strain and tilt records. These observational data have been used in rectangular fault models with uniform slip to analyze SSEs; however, the assumption of uniform slip precludes the possibility of mapping the slip distribution in detail. We report here an inversion method, based on the joint use of strain and tilt data and evaluated in terms of the Akaike’s Bayesian information criterion (ABIC), to estimate the slip distributions of short-term SSEs on the plate interface. Tests of this method yield slip distributions with smaller errors than are possible with the use of strain or tilt data alone. This method provides detailed spatial slip distributions of short-term SSEs including probability estimates, enabling improved monitoring of their locations and amounts of slip.


Introduction
Several types of slow earthquakes commonly occur along the Nankai Trough subduction zone ( Fig. 1) of the Philippine Sea plate southwest of Japan (Obara and Kato 2016). Among these, short-term slow slip events (SSEs) with durations of several days have been detected from records of borehole strainmeters (Kobayashi et al. 2006) and tiltmeters Hirose and Obara 2005). Strain and tilt data have been combined to represent these short-term SSEs on the basis of rectangular fault models with uniform slip (e.g., Itaba et al. 2012). In another approach, Nishimura et al. (2013) detected short-term SSEs from Global Navigation Satellite System (GNSS) data by fitting linear functions with and without offsets and estimated fault model parameters. Because these estimates assumed a single rectangular fault model with uniform slip and used a least-squares evaluation of error, the spatial distribution of slip may be biased.
Alternatively, Kano and Kato (2020) estimated the detailed cumulative slip distribution of short-term SSEs along the Nankai Trough subduction zone from 2004 to 2009 by stacking GNSS data. Hirose and Kimura (2020) used tilt-offset measurements to estimate slip distributions of short-term SSEs from January 2001 to March 2019 in the Shikoku region and gain insight into the slip budget along the strike direction of the Nankai Trough. Although these estimates of the slip distribution on the plate boundary surface yield more detailed slip locations than rectangular fault models, they employed GNSS or tilt data alone. In this study, we developed a joint inversion of strain and tilt data using the Akaike's Bayesian information criterion (ABIC; Akaike 1980) to estimate detailed slip distributions of short-term SSEs. We report the use of this method to calculate detailed spatial slip Open Access *Correspondence: t.tsuyuki@mri-jma.go.jp 1 Meteorological Research Institute, 1-1 Nagamine, Tsukuba, Ibaraki, Japan Full list of author information is available at the end of the article distributions of short-term SSEs and evaluate the error of the obtained slip distributions.

Inversion of slip distribution
We represent the observation vector as o = (o 1 , . . . , o N ) T , where N is the number of observational data components, in this case strain and tilt steps, and (. . . ) T denotes the transpose of the matrix. We seek to determine the slip distribution assuming that the changing trends in strain and tilt data are caused by non-uniform slip on the plate interface. For this purpose, the plate interface is divided into M grid points, and we assume a small subfault centered on each point that has a rectangular shape and uniform slip. where G is a Green's function coefficient matrix of dimensions N by 2M when a uniform unit slip is given to the small subfault at each grid point. The matrix components can be obtained from surface deformation due to a fault slip in a homogeneous elastic half-space, computed (1) o 1 o 2 = Gs + e 1 e 2 , by the method of Okada (1992). e 1 and e 2 are error vectors of tilt and strain data, respectively. It is assumed that e 1 and e 2 follow a normal distribution with 0 mean and variance-covariance matrix σ 2 1 E 1 and σ 2 2 E 2 , respectively. σ 2 1 and σ 2 2 are scale factors for the variance of e 1 and e 2 . No correlation is assumed between components, and the covariance elements in the variance-covariance matrix are set at 0. It is not appropriate to use the same scale factor for the variance of different observation data. Therefore, we introduce the hyper-parameter η 2 (= σ 2 2 /σ 2 1 ) to evaluate errors on different forms of observational data. This hyper-parameter η 2 expresses the relative weight of the strain and tilt datasets in the probability distribution (Funning et al. 2014). Then, the probability density function of the observation vector given a slip distribution s is represented as where matrix E is a function of η 2 = σ 2 2 /σ 2 1 , which is an N by N dimension matrix, and . . . represents the determinant of the matrix. We represent a slip vector s using a basis function in which the basis function is the boxcar function The slip distribution s is a smooth distribution represented by the 2M by 2M matrix L and 2M-dimension vector δ as (Fukuda and Johnson 2008) where L represents a discrete Laplacian and δ follows a normal distribution with 0 mean and the variancecovariance matrix ρ 2 I (where I is the unit matrix). The resulting slip distribution also follows a normal distribution with 0 mean and variance-covariance matrix L −1 ρ 2 L −1 T , and the probability density function of s is written as follows with R = L T L: (6) Ls = δ, According to Bayes' theorem, the probability density function of the slip vector for the observation vector o is obtained with a normalizing coefficient c to set the probability to 1 as where Following Yabuki and Matsu'ura (1992), to find the slip vector s and hyper-parameters σ 2 1 , η 2 , and α 2 to maximize this probability function, first, we hold fixed the three hyper-parameters. Then, the probability density is maximized when the slip vector s is With this s * and s(s * ) , we can write (7) p s|ρ 2 = 2πρ 2 −M �R� 1/2 exp − 1 2ρ 2 s T Rs .
(9)  From this equation, the slip distribution s is a normal distribution with a mean of s * and covariance σ 2 1 G T E −1 G + α 2 R −1 . Next, with the marginal likelihood of σ 2 1 , η 2 , and α 2 when a slip distribution is given We then use the ABIC (Akaike 1980), which is to be minimized based on the entropy maximization principle (Akaike 1977) To do this requires ∂P/∂σ 2 1 = 0 , which follows: Substituting this into (14), we obtain where C is a constant independent of α 2 and η 2 . As we cannot calculate α 2 and η 2 to minimize ABIC analytically, we determine optimal values α 2 and η 2 by the shuffled complex evolution (SCE-UA) method (Duan et al. 1992).
With these in hand, the optimal value of s is given as the optimal value of σ 2 1 is given as. and the covariance of the estimated errors of s is given as

Estimation of rectangular fault model parameters
For the evaluation of the slip-distribution model, we also modeled a short-term SSE as a rectangular fault with parameters where λ is the longitude, φ the latitude, and h the depth of the center of the fault; L is the length, W the width, S the strike, and δ the dip of the fault; and R is the rake and D the amount of slip. The amount of slip is uniform on the fault in this model. Observation data (o 1 , o 2 ) T are represented with residual r as where the Green's function coefficient matrix G is calculated by the same method used for the inversion of slip distribution. We estimate parameters of that minimize L2 norm First, the optimal position ( , φ ) of the fault and the amount of slip D are calculated by a grid search method with the other parameters held fixed, using the upper surface of the subducting plate as the position according to Itaba and Ando (2011). Since we fixed the length and the width of the fault 20 km in the grid search, strain and tilt are linearly related to the amount of slip. Second, the optimal value of parameter is calculated to minimize the sum of squares of Also shown are numbers of LFTs per hour and barometric pressure and precipitation data from Hamamatsu weather station. Linear trends have been removed from the original strain and tilt data. Data from the period in blue (0000 on 14 April to 0000 on 19 April) were used to estimate slip distribution in this study the residual between theoretical and observed values within a certain search range around the optimal fault position and fault parameter at that position by the SCE-UA method (Duan et al. 1992).

Strain and tilt data
We used strain observational data from volume strainmeters (Sacks et al. 1971) and multi-component strainmeters (Ishii et al. 2002) of the Japan Meteorological Agency (JMA), Shizuoka Prefecture, and the National Institute of Advanced Industrial Science and Technology (AIST). We used tilt data from the high-sensitivity accelerometers of the Hi-net (Okada et al. 2004;Obara et al. 2005)  ig. 8 Distribution of the standard deviation (1σ) of the estimated slip (Fig. 7) in the N170E direction (left panel) and N80E direction (right panel) at each subfault. Map region corresponds to Area A in Fig. 1. a, b Standard deviation of the estimated slip based on strain and tilt data. c, d Standard deviation of the estimated slip based on strain data alone. e, f Standard deviation of the estimated slip based on tilt data alone and geomagnetism were removed from the 1-h average strain and tilt values during the period from January 2015 to December 2016 as described below. Background linear trends during a time period before the SSEs were then removed from the corrected data.
Barometric pressure effects were corrected by first-order correlation with data from the barometer attached to the observation station or barometric data from the nearest JMA weather station. Tidal effects were removed from the pressure-corrected data using Baytap-G software (Tamura et al. 1991). Precipitation effects were corrected using a three-stage tank model . Because the multi-component strainmeters contain a magnetic sensor, the influence of geomagnetism was removed by the method of Miyaoka (2011). In addition, observed (in-situ) strain data from borehole strainmeters may be amplified or attenuated in response to crustal strain (Gladwin and Hart 1985). This effect was removed by a calibration based on theoretical tides (Matsumoto et al. 2001;Hirose et al. 2019) for multi-component strainmeters and one based on theoretical teleseismic waveforms (e.g., Nakanishi and Kasahara 1990) for volume strainmeters.

Subfaults parameters
The set of subfaults used for the basis functions introduced in the Method section is shown in Fig. 2. We  Fig. 1. a Strain data and tilt data from observations (blue) and calculated from the estimated slip distribution shown in Fig. 7a (white). b Strain data from observations (blue) and calculated from the estimated slip distribution shown in Fig. 7b (white). c Tilt data from observations (blue) and calculated from the estimated slip distribution shown in Fig. 7c (white) approximated the area around the Nankai Trough by rectangular coordinates with the Y-axis oriented N55°W, the direction of subduction of the Philippine Sea plate (Heki and Miyazaki 2001). Grid points were set every 10 km in the X-axis (N35°E) direction and 10 km intervals along the plate interface. X and Y axes are shown in Fig. 2. We defined a small 10 km square subfault (rectangular in Fig. 2) centered on each grid points with specified strike and dip according to the configuration of the upper surface of the subducting Philippine Sea plate (Hirose et al. 2008) and used in subsequent analyses. As we represented a slip vector using the boxcar basis function (Eqs. 4, 5), we supposed a uniform slip at each subfaults and calculated the Green's function. It should be mentioned that subfaults have overlaps and gaps on the plate interface in discretization, which can cause errors when calculating the whole slip amount in the region.

Verification of the method using a synthetic dataset
We checked the validity of our inversion method using a synthetic dataset for the set of subfaults in region A in Fig. 2. The subfaults were assigned the slip distribution shown in Fig. 3a,b, in which the slip amount varies from point to point (every 10 km). We calculated synthetic linear strain at 14 measurement points and tilt at 10 measurement points in a homogeneous elastic half-space with a Poisson's ratio of 0.25 using the formula of Okada (1992), yielding the results shown in Fig. 4. Strain changes were converted to the principal strain in this figure. We did not consider the measurement errors. We estimated a slip distribution (Fig. 3c,d) and rectangular fault model parameters separately with this synthetic data set. Since we did not consider the measurement errors, only modeling errors were contained (Yabuki and Matsu'ura 1992) in the errors of Eq. (1). Therefore, in the variancecovariance matrix of observation error, the covariance components were not considered and the variance components were set to 1.0.
Optimal hyper-parameters from this exercise are listed in Table 1. We introduced hyper-parameter η 2 to represent the weight of the strain and tilt dataset. In this test, since the observation data used for the inversion do not include any error, this parameter representing the weights of the strain and tilt data should be 1.0. Therefore, it is reasonable to obtain 0.92 as this value. The slip distribution retrieved by the inversion ( Fig. 3c, d) nearly matched the given slip distribution (Fig. 3a, b). Distribution of the standard deviation of the estimated slip in the two orthogonal directions at each subfault is shown in Fig. 3e, f. The standard deviation was used as a reference to evaluate the estimation error of slip obtained in each subfault. It can be seen that the standard deviation is relatively small in the region where many observation points are distributed. We calculated the total seismic moment using only subfaults with slips larger than the 1-σ estimation error in this study. The calculated total seismic moment of 1.45 × 10 18 N m (M w = 6.04) was almost equal to the given seismic moment of 1.41 × 10 18 N m (M w = 6.03). A rigidity of 4.0 × 10 10 N m -2 was assumed in this study. Strain and tilt change calculated from the given and the retrieved slip distribution are similar (Fig. 4). We also estimated a rectangular fault model with the same data set. The estimation based on the rectangular fault model (rectangle in Fig. 3d) imperfectly reproduced the given slip and the eastern side of the given slip area is not included in the rectangular fault model. Estimated seismic moment was 1.33 × 10 18 N m (M w = 6.02). Slip on the fault was assumed to be uniform in rectangular fault model, and therefore, the estimated fault model, while including all points with given large slip, excluded some with given small slip. When calculating the average slip amount for a long period, the estimation of the slip distribution at the grids is more appropriate than the rectangular fault.

Results and discussion
In this section, we analyzed some known short-term SSEs as examples to test our inversion method.

Example 1
In the first example, we tested whether the simultaneous use of strain and tilt data (joint use) improves the estimation error for a short-term SSE beneath region A in Fig. 2. Figure 5 shows the distribution of low-frequency tremor events (LFTs) in the Tokai area during 14-19 April 2015. We plotted the LFTs taken from the NIED catalog (Maeda and Obara 2009;Obara et al. 2010) in this study. During that time, the strain and tilt change appeared (Fig. 6). These phenomena may correspond to a short-term SSE in this region (e.g., . We analyzed the slip distributions resulting from the use of strain or tilt data alone, and analyzed the slip distribution by our inversion method using strain and tilt data. The observational components to be used were selected manually. Observed values are calculated using the difference between the corrected data at the beginning and the end of the period in this study (e.g., Itaba and Ando 2011). The uncertainties of the data are represented by the variance from 11 March to 14 April 2015. (Additional file 1: Table S1). We considered that the error contained in the data was mainly observation error, and that the observation component with larger variation in the original data had larger error. Therefore, in the variance-covariance matrix, the diagonal elements were set to the ratio of each component's variance to their maximum variance for the  strain and tilt data respectively. The optimal hyperparameters are listed in Table 1. Estimated slip distribution, distribution of the standard deviation of the slip in the two orthogonal directions at each subfault, and the comparison between observed strain and tilt data and calculated strain and tilt data are shown in Figs. 7,8,9,respectively. In the comparison between observed strain and calculated strain, linear strain components of each station used for inversion were converted to the principal strain. Seismic moment and the standard deviation and amount of the slip in the two orthogonal directions at the maximum slip subfault are listed in Table 2. The standard deviation (1σ estimation error) is larger when using stain data or tilt data alone than by joint use (Fig. 8). The ratio of the standard deviation to the estimated slip amount in the two orthogonal directions at the maximum slip subfault are 77% and 51% when using strain data alone, while 16% and 34% for joint use (Table 2). We could not estimate slip larger than the standard deviation when using tilt data alone (Fig. 7c). Also, when using strain data alone, slips in the same direction as the plate convergence are shown in the eastern side of the grids which are not considered to be due to the short-term SSE (Fig. 7b). The addition of tilt data therefore yielded a more plausible slip distribution with better estimation error.
We denoted the relative weight of strain and tilt data as a value of hyper-parameter η 2 . The value of η 2 is 0.38 (Table 1) in this example, signifying that the contribution of the strain data was relatively large. Because the average variance ratio of strain and tilt data used in this study was about 0.5, the ratio obtained by inversion was slightly small.

Example 2
In the second example, we tested whether the slip distribution can be tracked as a short-term SSE migrates beneath region A in Fig. 2. During 29 December 2015-6 January 2016, LFTs shifted toward the northeast in the Kii Peninsula and the Tokai area (Fig. 10). At the same time, both strain and tilt data change appeared (Fig. 11). These phenomena may correspond to a short-term SSE in this region (e.g., . We divided the data into three periods (Fig. 11) and analyzed slip distributions by our inversion method using strain and tilt data during each period. Observed values and their observation errors are listed in Additional file 1: Table S2. In the variance-covariance matrix, the diagonal elements were set to the ratio of each component's variance to their maximum variance for the strain and tilt data, respectively, from the preceding period (15-29 December 2015). The optimal hyper-parameters are listed in Table 1. Estimated slip distribution, distribution of the standard deviation of the slip in the two orthogonal directions at each subfault, and the comparison between observed strain and tilt data and calculated strain and tilt data are shown in Figs. 12, 13, 14, respectively. We also estimated rectangular fault model parameters (rectangle in Fig. 12) with the same dataset. Seismic moment and the standard deviation and amount of the slip in the two orthogonal directions at the maximum slip subfault are listed in Table 2. The slip distribution obtained by inversion follows the propagation of the LFTs and they probably indicate the migration of the short-term SSE. Estimated distribution of slips larger than 1σ estimation error is almost same as the distribution of LFTs, whereas the slip distribution of the first period is slightly wide. This may be because there is no observation point on the west side of the slip. The locations of the rectangular faults were generally consistent with the LFTs, but it was estimated to be a bit wider for the third period.

Example 3
In the third example, we analyzed a short-term SSE that occurred in the Bungo Channel during 23-27 October, 2016 (e.g., Ochi et al. 2017;Kimura 2017) where LFTs were also observed (Fig. 15). Both the strain and tilt data appeared to respond to this event (Fig. 16). We analyzed the slip distribution using strain and tilt data during the period shown in color on Fig. 16 (2100 on 23 October to 0000 on 27 October) and the subfaults of region B in Fig. 2. We also analyzed rectangular fault model parameters using the same data set. Strainmeters are not as densely distributed in this region as in the Tokai region (Fig. 1). Therefore, we also analyzed the slip distributions resulting from the use of strain or tilt data alone and compare the results with the one obtained with joint use. Observed values and their observation errors are listed in Additional file 1: Table S3. The diagonal elements in the variance-covariance matrix were set to the ratio of each component's variance to their maximum variance for the strain and tilt data, respectively, from 10 to 23 October 2016. The optimal  Fig. 1. a Strain and tilt data for the first period (29 December 2015 to 1 January 2016) from observations (blue) and calculated from the estimated slip distribution shown in Fig. 12a (white). b Strain and tilt data for the second period (1 January to 1200 on 5 January) from observations and calculated from the estimated slip distribution shown in Fig. 12b. c Strain and tilt data for the third period (1200 on 5 January to 7 January) from observations and calculated from the estimated slip distribution shown in Fig. 12c   hyper-parameters are listed in Table 1. Estimated slip distribution, distribution of the standard deviation of the slip in the two orthogonal directions at each subfault, and the comparison between observed strain and tilt data and calculated strain and tilt data are shown in Figs. 17, 18, 19, respectively. Seismic moment and the standard deviation and amount of the slip in the two orthogonal directions at the maximum slip subfault are listed in Table 2. Although the standard deviation of the estimated slip is small (Fig. 18) in the result using strain data alone, small slips are widely distributed (Fig. 17b) and there is almost no slip in the distribution of LFTs. When using tilt data alone, we could not estimate slip larger than the standard deviation (Fig. 17c). On the other hand, the standard deviation of each subfault slip estimated by joint use of strain and tilt data is smaller than when using tilt data alone. Region of subfaults where the estimated slip is larger than 1σ estimation error is almost the same as the distribution of LFTs (Fig. 17a). In this example as well as Example 1, we could obtain more likely slip distribution by joint use of strain and tilt data.

Conclusion
Short-term SSEs along the subduction zone of the Nankai Trough are typically detected by interpreting strain and tilt data. We formulated an inversion method to estimate the slip distribution of these SSEs by the joint use of strain and tilt data based on ABIC. The effectiveness of this method was demonstrated by a test using a synthetic dataset. When we applied this method to analyze the slip distributions of known short-term SSEs, we showed that estimated slip distribution follows the migration of the short-term SSE. We also showed that a plausible slip distribution with less estimation error can be obtained by the joint use of strain and tilt data than when each data is used alone. By stochastically evaluating the slip distributions obtained by our method, we showed that it is possible to map variations in the detailed slip distribution on the plate boundary, an improvement over results using rectangular fault models with uniform slip.  Fig. 1. a Strain data and tilt data from observations (blue) and calculated from the estimated slip distribution shown in Fig. 17a (white). b Strain data from observations (blue) and calculated from the estimated slip shown in Fig. 17b (white). c Tilt data from observations (blue) and calculated from the estimated slip distribution shown in Fig. 17c (white)