On the use of tsunami-source data for high-resolution fault imaging of offshore earthquakes

The source imaging for offshore earthquakes using terrestrial geodetic data has a limited estimation performance due to the low data resolution. One approach to overcome this limitation is the use of seafloor geodetic data. In this study, we focus on tsunami-source data, which is the spatial distribution of vertical crustal displacements above the source area and can be derived from tsunami waveform records. We evaluate how the use of this spatial seafloor geodetic data improves the estimation of a rectangular fault model. Here, the fault model of the 2016 off-Fukushima earthquake in Japan, which was a shallow intraplate earthquake (Mw 7.0), was estimated by three inversions: terrestrial Global Navigation Satellite System (GNSS) data only, tsunami-source data only, and a combination of the GNSS data and tsunami-source data. A Bayesian inversion approach was used to understand the distribution of the estimated fault parameters and their relationship. The results indicated that the terrestrial GNSS data have a low resolution for the analysis of the offshore earthquake, which resulted in a biased solution with large uncertainty. Conversely, the use of tsunami-source data significantly improved the resolution and reliability of source imaging and reduced the dependency among fault parameters. These results suggested that the high-spatial-resolution information of tsunami source is a powerful tool in source imaging of offshore shallow earthquakes. Moreover, the combined use of the two different geodetic data leads to a more robust estimation of fault parameters. We believe that the use of tsunami-source data is useful, not only for the post hoc source analysis, but also for estimating an earthquake rupture area just after a large earthquake, where GNSS data are currently used.


Background
Global Navigation Satellite System (GNSS) has been deployed over a wide land area, and the GNSS-based static crustal deformation data have played an important role in estimating fault rupture for large earthquakes (e.g., Ozawa et al. 2011;Crowell et al. 2012).However, source imaging using only terrestrial geodetic data has limited resolution for earthquakes occurring in offshore areas (Romano et al. 2010;Kubo and Kakehi 2013;Yoshioka and Matsuoka 2013;Williamson and Newman 2018).The use of seafloor geodetic data has been used to address this issue; previous studies have demonstrated that the use of oceanic geodetic data improves the resolution of the source imaging (Kubo and Kakehi 2013;Yoshioka and Matsuoka 2013;Yokota et al. 2016).There are several types of seafloor geodetic data, such as the crustal movement obtained by a combination of GNSS and acoustic ranging (GNSS-A) technique (e.g., Spiess et al. 1998;Sato et al. 2011;Yokota et al. 2018) and vertical displacements directly observed by ocean-bottom pressure gauges (OBPG) (e.g., Mikada et al. 2006;Ito et al. 2011).In this study, we focus on the seafloor geodetic data inferred from tsunamis.The tsunami source, which is the spatial distribution of vertical crustal displacements above the source area of a large earthquake, has been estimated from tsunami waveform records for large earthquakes (Aida 1972;Baba et al. 2005;Tsushima et al. 2009Tsushima et al. , 2011;;Wei et al. 2013;Dettmer et al. 2016;Kubota et al. 2018Kubota et al. , 2021)).Most observations of GNSS-A and OBPG are sparsely distributed far from an earthquake source area; therefore, available seafloor geodetic data are often limited, and its use may not improve the resolution of fault imaging as much as expected.On the other hand, the tsunami source provides rich information on the earthquake fault rupture, especially for shallow earthquakes, because it is a spatial distribution of seafloor vertical displacement above the source area.The use of tsunami-source data is useful not only for the post hoc analysis including this study, but also for the analysis just after large earthquakes.
The purpose of this study is to evaluate how the additional use of tsunami-source data improves the estimation of a rectangular fault for offshore earthquakes, compared with using only the terrestrial GNSS data alone.Although the estimation of a heterogeneous faultslip distribution based on geodetic data is often done, the estimation of a rectangular fault is essential for understanding of the entire source image including the extent of fault rupture area because this estimation can provide a stable solution without assuming the fault plane.Moreover, the estimated result just after the occurrence of a large earthquake is useful for the forward simulation of tsunami early warning.Although the tsunami waveform data can be directly used in an inversion to estimate the fault model (e.g., Satake 1989;Gusman et al. 2012;Kubota et al. 2021Kubota et al. , 2022)), we adopt the use of the tsunami source.This approach improves the solution stability in the quality assurance of tsunami source estimation and the consistency of the tsunami-driven seafloor geodetic data with the GNSS data (crustal deformation data).Our approach is similar to the two-step tsunami source inversion method of Gusman et al. (2018a, b) in which tsunami waveforms are firstly inverted for the initial seasurface height distribution and the estimated tsunami source is then inverted for the fault slip distribution.For the inversion method, we adopt a Bayesian approach to evaluate the uncertainty of estimated fault parameters and the dependence among the fault parameters.
For this purpose, we estimate the rectangular fault model of the 2016 off-Fukushima earthquake by the joint inversion using both terrestrial GNSS data and tsunamisource data and the inversions using only one of the two data, and compare these results.This earthquake occurred off the Fukushima Prefecture, Tohoku, in eastern Japan at 05:59 on November 22, 2016, (JST; 20:59 on November 21, UTC).The Japan Meteorological Agency (JMA) magnitude (M JMA ) of this event was 7.4, and the moment magnitude (M w ) provided by moment tensor analysis of F-net of National Research Institute for Earth Science and Disaster Resilience (NIED) (Fukuyama et al. 1998) was 7.0.The hypocentral location, focal mechanism, and spatial distribution of aftershocks indicated that this earthquake was a shallow normal-faulting crustal earthquake within the overriding plate (Fig. 1).The consequent tsunami caused by this earthquake was observed along the coast of Tohoku and Kanto regions, and its maximum amplitude was 1.4 m at the Sendai port in Miyagi Prefecture (Japan Meteorological Agency 2016).

Data and methods
Herein, we conducted three different inversions, namely joint inversion using both terrestrial GNSS data and tsunami-source data (joint inversion), inversion using only GNSS data (GNSS inversion), and inversion using only tsunami-source data (tsunami-source inversion).Hereafter, we comprehensively explain joint inversion.GNSS inversion and tsunami-source inversion follow almost the same procedure, except that a hyperparameter weight between different data is not used and the error scale factor of each observation equation is pre-set to the same value as the joint inversion.
For the terrestrial GNSS data, we used two horizontal components of static displacements at 75 stations of the GNSS Earth Observation Network System (GEONET) of the Geospatial Information Authority of Japan (GSI) (Yamagiwa et al. 2006;Sagiya 2004) (Fig. 1).These were obtained by calculating the difference between the preseismic and postseismic positions of daily coordinates of GEONET (F3 solution) (Nakagawa et al. 2009), which were acquired by averaging the positions from November 14 to November 20 (JST) and from November 22 to November 26 (JST), respectively.Static horizontal displacements toward the east or northeast direction with a maximum of approximately 0.05 m were found at the GNSS stations along the coast of Fukushima Prefecture close to the hypocenter (Fig. 1).
For the tsunami-source data, we used the spatial distribution of seafloor vertical displacements above the source region (Fig. 1) estimated by Kubota et al. (2021).Based on the tsunami waveforms of the OBPG records at stations of NIED S-net (Kanazawa et al. 2016;Mochizuki et al. 2016;Uehira et al. 2016;Aoi et al. 2020), Kubota et al. (2021) inverted the spatial distribution of initial sea-surface height (tsunami source) with the horizontal spatial interval of 2 km, in an area of 50 km × 50 km with the imposition of the smoothing constraint on its distribution.A subsidence with a horizontal extent of ~ 40 km × ~ 20 km along the northeast-southwest direction, having a peak of ~ 2 m, was found (Fig. 1).
We estimated a rectangular fault model with a homogeneous slip following Kawamoto et al. (2016Kawamoto et al. ( , 2017) ) and Ohno et al. (2021).Unknown fault parameters were latitude, longitude, and depth at the top center of the fault, strike angle, dip angle, rake angle, fault length, fault width, and fault slip.To acquire the fault parameters of a single rectangular fault model through the Bayesian framework, we used the Markov Chain Monte Carlo (MCMC) method following Ohno et al. (2021).We assumed the observation equation wherein fault parameters m ( M-dimensional vector) are related to static dis- placements at receivers d ( N-dimensional vector) via Green's functions G ( N × M matrix): where N and M are the numbers of data and model parameters, respectively, and ε is an N-dimensional error vector.For simplicity, we assumed that the errors ε fol- low a Gaussian distribution and are independent of each other, ε ∼ N (0, σ 2 I N ) , where I N is an N-dimensional identity matrix and σ 2 is a scale factor of the errors ε .The probability density function (PDF) of the observation equation was as follows: where �a� 2 = a T a .σ 2 is also called hyperparameter.As σ 2 must be positive, σ 2 is transformed into 10 Σ and Σ is sampled instead of σ 2 (Kubo et al. 2016):   If two data, d 1 and d 2 , which are independent of each other, are used, the PDF of the observation equations is as follows: where N 1 and N 2 are the number of data for d 1 and d 2 , respectively, G 1 and G 2 are the Green's functions corre- sponding to d 1 and d 2 , respectively, and Σ 1 and Σ 2 are the scale factors of the errors for d 1 and d 2 , respectively.Connecting the observation equation and the prior information using Bayes' theorem (e.g., Tarantola 2005), the conditional PDF of the unknown parameters was related to the product of the PDF expressing prior information on the unknown parameters (prior PDF) by the PDF of the data given the unknown parameters (PDF of the observation equation): Here, unknown parameters are the fault parameters of rectangular fault ( m ) and scale factors of errors for terrestrial and seafloor geodetic data ( Σ 1 and Σ 2 ).The unknown parameters were estimated as follows: (1) for several values of the relative weight between the two datasets ( A = Σ 2 /Σ 1 ), the posterior PDF of m and Σ 1 was ensembled using the basic Metropolis-Hasting (M-H) method (Metropolis et al. 1953;Hastings 1970) (5) determined to ensure that the data fit for each dataset was satisfactory; (3) the posterior PDF of m and Σ 1 for the fixed value of A was ensembled.The prior PDFs of m and Σ 1 were assumed to follow a uniform distribution (Table 1).The initial parameters of the fault parameters were set by following the moment tensor solution of NIED F-net and the scaling relationships among fault parameters (Table 1).The width of the proposal probability density, which controls the efficiency of the M-H method, was set by trial and error following Gelman et al. (1996).Although the results of the Bayesian inversion are represented as the posterior PDF of unknown parameters, the maximum a posteriori (MAP) solution was obtained as the optimal fault model as follows: (1) the posterior PDF of m and Σ 1 was ensembled by the M-H method; (2) the median of the marginal posterior distribution for Σ 1 was determined as its optimal value; (3) given the optimal value of Σ 1 , the posterior PDF of m was ensembled by the M-H method; (4) the combination of the model parameters with the maximum posterior probability was adopted as the optimal fault model.For Green's functions, we calculated the theoretical static displacements caused by a unit slip on the rectangular fault, assuming a homogeneous elastic half-space (Okada 1992).Poisson's ratio and rigidity were assumed to be 0.25 and 30 GPa, respectively.An MCMC

Initial values Prior PDFs
Latitude (°) 37.3547 a U(−180, 180) sampling with 10 6 steps was conducted, the first 2.0 × 10 5 steps were discarded as burn-in, and every 10th model visited in the last 8.0 × 10 5 steps was taken to obtain 8.0 × 10 4 ensembles as the posterior inference.
The computation time required for the MCMC sampling with 10 6 steps is ~ 110 min when using a single processor of Xeon E5-2630 v3 (2.4 GHz).

Results
The results of the joint, GNSS, and tsunami-source inversions were compared.First, we focused on their optimal fault models.Table 2 lists the optimal fault models in the three inversions, which are also shown in Fig. 2. Although, both inversions suggested a fault model with the southeast-dipping normal-faulting fault located southwest of the epicenter, the three models showed some differences.Compared with the joint and tsunamisource inversion models, the GNSS inversion model had a longer fault length and a much shorter fault width.Moreover, the estimated slip value of the GNSS inversion model was much larger than that of the other models, and its strike was rotated counterclockwise and the dip was much gentle.Although the other two inversion models were relatively similar, the estimated fault of the tsunami-source inversion model was rotated clockwise compared with that of the joint inversion.Figure 2 indicates that the fault of the GNSS inversion model did not overlap with the spatial distribution of aftershocks by JMA.The strike angle of the GNSS inversion model was gentler than that of the south-east-dipping nodal plane of the F-net moment tensor solution (Fig. 2 and Table 1).
The joint and tsunami-source inversion models were consistent with the spatial aftershock distribution and the F-net moment tensor solution.
The data fittings of the optimal models for terrestrial and seafloor geodetic data are shown in Figs. 2  and 3, respectively, and their variance reduction values are listed in Table 2.Although the GNSS and tsunamisource inversion models showed a better data fitting than the joint inversion model for each used data, the joint inversion model also reproduced both observations sufficiently.The seafloor geodetic deformation caused by the GNSS inversion model was localized, and their spatial distribution was inconsistent with the smooth seafloor geodetic distribution of Kubota et al. (2021), as shown in Fig. 3c, f.This indicates that the seafloor geodetic distribution of Kubota et al. (2021) implies the deeper fault with a steeper dip and a smaller slip than the GNSS inversion model.The synthetic terrestrial static displacements of the tsunami-source inversion model were not consistent with the observations (Fig. 2c).In particular, the direction of the displacement vector in the area denoted by a broken circle in Fig. 2c was different between the observations and synthetics, which indicates that the strike angle of the tsunami-source inversion model was rotated too far clockwise.
Next, we focused on the posterior PDF for fault parameters in the three inversions.Figure 4 shows the relationship among the estimated fault parameters in their ensemble of the joint inversion.The histograms of the fault parameters suggested that the marginal posterior PDF of each fault parameter had a simple unimodal distribution.The heatmaps among fault parameters indicated several dependencies among fault parameters.For example, longitude, latitude, strike angle, and rake angle are correlated.Because the seafloor subsidence zone extended along the northeast-southwest direction (Fig. 3a), the estimated horizontal location of the rectangular fault varied along the northeast-southwest direction.The variations in the fault position could cause variations in the strike angle and rake angle, which depend on each other due to the tradeoff involved in decomposing the horizontal slip vector into strike and rake angles.Moreover, the dip angle and fault width were negatively correlated.The horizontal location of the bottom edge of the rectangular fault could be constrained by the spatial distribution of seafloor crustal deformation; however, the dip angle and fault width inferred from the horizontal information of the bottom edge were not unique, resulting in a negative correlation between the two fault parameters.The fault slip was negatively correlated with the fault length and width.The top depth of the fault was correlated with the slip, probably because of their tradeoff for the effect on surface crustal deformation.The relationship between the variance reduction value of each data and strike angle indicates the tradeoff of their data fittings with strike angle; the clockwise rotation of fault strike can better explain the terrestrial GNSS data (especially the area denoted by a broken circle in Fig. 2), whereas the counterclockwise rotation of fault strike can better explain the seafloor geodetic data of Kubota et al. (2021).
Figure 5 shows the comparison of the relationship of the fault parameters among the joint, GNSS, and tsunami-source inversion models.The fault-parameter relationship of GNSS and tsunami-source inversion models are shown in Additional file 1: Figure S1, S2, respectively.The variation in the estimated fault parameters in the GNSS inversion was greater than that in the joint inversion.Moreover, the relationship among the fault parameters in the GNSS inversion had a significantly complex dependency on each other, indicating poor resolution in analyzing offshore earthquakes using only terrestrial data.In the GNSS inversion, the latitude and longitude were significantly negatively correlated, and this direction of the location uncertainty roughly corresponded to the direction of the horizontal slip vector.This implies that the spatial location of the rectangular fault was not resolved along the northwest-southeast direction in the GNSS inversion.This was because the GNSS stations were unevenly distributed for offshore earthquakes, which consequently led to the low spatial resolution along the slip-vector direction.This correlation trend between latitude and longitude was not found in the joint and tsunami-source inversion models because the seafloor geodetic data suppressed the positional uncertainty along the northwest-southeast direction.We also found that the dip angle was positively correlated with the fault slip in the GNSS inversion.In general, the horizontal static displacements are useful to constrain the horizontal components of fault slip; however, it is difficult to decompose the horizontal components of fault slip into the dip angle and the net slip on the fault, especially in the case where the horizontal geodetic observations are far located from the source area.This is found as the dependency between dip angle and fault slip in the GNSS inversion model.One approach to this problem is to use the vertical component of geodetic data, such as the tsunami source; therefore, the dependency between the dip angle and fault slip is not found in the joint and tsunami-source inversion models.The trend of fault-parameter dependence in the tsunami-source inversion was similar to that in the joint inversion.The variation of the fault parameters in the joint inversion was slightly smaller than that in the tsunami-source inversion.
Figure 6 shows the spatial variation of the estimated fault in the three inversions.The solutions in the joint and tsunami-source inversions almost converged, whereas the solutions in the GNSS inversion varied widely.The marginal distribution of the fault location in the GNSS inversion was inconsistent with the spatial distribution of the aftershocks.

Advantage of using tsunami-source data for fault imaging
This study demonstrated that the posterior PDFs of some fault parameters in the GNSS-only inversion were biased.Their fault models had shorter widths and longer lengths with a shallower top depth than those in the other inversion models.Furthermore, the GNSS-only inversion's model had a counterclockwise-rotated strike angle, a more gentle dip angle, and a much larger slip compared with those of the joint and tsunami-source inversion.Moreover, the spatial location of the fault model shifted northwestward.The features of the fault model derived from only terrestrial GNSS data are inconsistent with other source information, such as the spatial distribution of aftershocks, the source mechanism solution, and source-process models estimated in previous studies (Gusman et al. 2017;Kubota et al. 2021).This indicates the difficulty in resolving the source information for offshore earthquakes based only on the terrestrial geodetic observations that are one-sidedly distributed far from the source region.This difficulty is also reflected in the significantly large uncertainty of the estimated fault parameters in the GNSS inversion.The posterior PDFs in the GNSS inversion suggest that the fault location was not resolved along the offshore direction from the terrestrial geodetic observations.This can occur with a one-sided distribution of observation stations located far from the source region; i.e., when the source of an offshore earthquake is analyzed using terrestrial geodetic data alone.The posterior PDFs in the GNSS inversion also suggest a strong correlation between dip angle and slip caused by the difficulty of decomposing dip angle and fault slip for an offshore earthquake based on horizontal components of terrestrial static displacements.This difficulty is expected to be pronounced for earthquakes with a large amount of dip-slip and a non-steep dip angle, which corresponds to the analysis of the 2016 off-Fukushima earthquake.The similar results would be true for other terrestrial geodetic data such as Interferometric Synthetic Aperture Radar (InSAR) data.
Fixing some fault parameters with other information is expected to improve the accuracy and uniqueness of the solution in analyses using only GNSS.For example, a centroid moment tensor (CMT) solution can provide information on fault mechanism (strike, dip, and rake angles) and the information is useful for obtaining reliable and robust solutions.However, not all tradeoffs among fault parameters can be eliminated by using the CMT solution, and the uncertainties of fault parameters, such as the spatial location of fault, its geometry, and fault slip would remain.Moreover, the CMT solution itself also has uncertainty.For example, there is a tradeoff between dip angle and seismic moment in the CMT solution for shallow earthquakes (Kanamori and Given 1981).
The results also suggested that the use of tsunamisource data could significantly improve the resolution and reliability of the source analysis for offshore earthquakes.We considered that the tsunami-source data are similar to the InSAR observation data as both data have a high spatial resolution for crustal movement in wide regional coverage.High-spatial-resolution geodetic data, such as InSAR data, can detect fine-scale crustal deformation and assist in comprehensively understanding earthquake fault ruptures (e.g., Massonnet et al. 1993;Biggs and Wright 2020).The tsunami-source data are useful for source-process analysis of offshore earthquakes, especially shallow intraplate earthquakes and interplate earthquakes with large slips near the trench, which are more likely to excite tsunamis.
The fault model in the tsunami-source inversion is close to that in the joint inversion; however, they differ in the strike angle, where the strike angle in the tsunamisource inversion is rotated too clockwise to reproduce the GNSS data.Although the tsunami sources are densely distributed, they are the vertical crustal deformation information.The additional use of other horizontal geodetic information is expected to lead to a more robust solution, which is supported by the relatively small variation of fault parameters in the joint inversion.

Application of a Bayesian approach to fault imaging
Previous studies that assessed the rectangular fault model have barely discussed the dependencies among parameters.Ohno et al. (2021) indicated that shallow offshore earthquakes result in multiple peaks in the marginal distributions of fault parameters and the dependency between strike and rake angles.Using the Bayesian approach, this study highlighted the complex relationship among estimated fault parameters in the source analysis for offshore earthquakes using only terrestrial geodetic data.Specifically, our results indicated that when badcondition inverse problems occur, investigating the entire distribution of solutions is necessary, rather than estimating only a single optimized solution.We also demonstrated that the use of seafloor geodetic data significantly improves the inverse problem and drastically reduces the uncertainty in the estimated fault parameters.

Potential of using tsunami-source data in rapid fault estimation
Following the recent development of GNSS networks that can provide real-time data on crustal deformation, the rapid estimation methods of fault parameters after large earthquakes using such real-time geodetic data have been developed (e.g., Ohta et al. 2012Ohta et al. , 2016;;Melgar et al. 2012;Crowell et al. 2012Crowell et al. , 2016;;Grapenthin et al. 2014;Kawamoto et al. 2016Kawamoto et al. , 2017;;Ohno et al. 2021).The accuracy of its estimation is important for disaster response; however, the resolution of GNSS data is limited for offshore earthquakes, as shown in this study.Our results indicate that the additional use of tsunami-source data is a good approach for this problem.The recent development of the offshore observation networks of OBPG such as S-net in Japan enables the quasi-real-time estimation of the tsunami source using real-time OBPG data (Titov et al. 2005;Baba et al. 2005;Tsushima et al. 2009Tsushima et al. , 2011;;Wei et al. 2013;Dettmer et al. 2016;Kubota et al. 2018Kubota et al. , 2021;;Suzuki et al. 2020).From the perspective of the quasi-real-time analysis, the tsunami-source data are superior to other seafloor geodetic data such as GNSS-A data and static displacement data observed by OBPG.Acquiring the GNSS-A records immediately after a large earthquake is currently challenging.Although the realtime pressure data recorded by OBPG includes information on static displacements due to an earthquake, real-time extraction of static displacement components from the OBPG data is sometimes difficult.For example, the S-net pressure waveform during earthquakes sometimes includes system-induced offset noise (Kubota et al. 2021) and it causes the difficulty to accurately separate the static displacement and the offset noises in real time.Moreover, the spatial resolution of tsunami-source data is superior to that of GNSS-A data and OBPG offset data, of which observation points are sparsely distributed.
Because of the tsunami propagation, the acquisition time of the tsunami source is expected to be several tens of minutes after an earthquake occurs (Kubota et al. 2018;Suzuki et al. 2020) and depends on the spatial location of the earthquake and the distribution of tsunami observation stations.Considering the acquisition timing of tsunami-source data, the fault information inferred from tsunami-source data would be useful not for the early warning of earthquakes and tsunamis, but for the identification of the rupture area.In the Nankai Trough subduction zone, great earthquakes (magnitude > 8) have occurred repeatedly at intervals of 100-200 years, and these great earthquakes have often occurred as pairs of M ~ 8 events occurring successively within an interval of 2 years (Ishibashi 2004).In such subduction zones with a history where a large earthquake has followed the most recent large earthquake, understanding which areas were ruptured in the most recent event is essential for assessing the probability of future fault rupture in the vicinity.
In this study, we focus on the estimation of a rectangular fault plane with a homogeneous slip.This estimation is expected to be useful for M7 earthquakes such as the 2016 off-Fukushima earthquake.On the other hand, it remains debatable whether this estimation approach is valid for huge earthquakes of M8 ~ 9 because their fault slip distributions were likely to be complex and may not be appropriated by a rectangular fault plane model with a homogeneous slip.Although the tsunami-source data are available even for huge earthquakes and would be useful for understanding the fault rupture process of huge earthquakes (Saito et al. 2011), it is necessary to validate this estimation approach through synthetic tests that simulate huge earthquakes.Another promising approach is to estimate the fault slip distribution on the pre-assumed fault that corresponds to a plate boundary (Kawamoto et al. 2017).

Conclusions
In this study, we evaluated how the use of tsunamisource data improves the estimation of rectangular fault model source imaging through the analysis for the 2016 off-Fukushima earthquake with a Bayesian inversion approach.Our results demonstrated that the terrestrial GNSS data have a low resolution for the analysis of the offshore earthquake, which resulted in a biased solution with large uncertainty.Conversely, the additional use of tsunami-source data significantly improved the resolution and reliability of source imaging and reduced the dependency among fault parameters.These results suggested that the high-spatial-resolution information of the tsunami source is useful in source imaging offshore shallow earthquakes.Moreover, the combined use of the two different geodetic data leads to a more robust estimation of fault parameters.We believe that the use of tsunamisource data is useful, not only for the post hoc source analysis, but also for the estimation of an earthquake rupture area just after a large earthquake, which is beneficial for the identification of the fault rupture area in the recent large earthquake.Although the existing systems for the rapid estimation of earthquake fault have used only GNSS data, we propose that the additional use of tsunami-source data can increase the analysis resolution for offshore earthquakes.Because there is a limitation in using data obtained from a single observation network, developing a strategic analysis that can successfully combine data obtained from each observation network is beneficial.

Fig. 1
Fig. 1 Map of the study area.The star indicates the epicenter of the 2016 off-Fukushima earthquake.The focal mechanism refers to the F-net moment tensor solution.The triangles indicate the GEONET stations used in this study.The arrows indicate the observed horizontal static displacements at the GEONET stations.The color distribution denotes the spatial distributions of uplift and subsidence estimated by Kubota et al. (2021).The circles represent the hypocenters of the aftershocks (M ≥ 3) within one day after the mainshock

a
Source location of the F-net moment tensor solution b Values of the south-east-dipping nodal plane of the F-net moment tensor solution c Values obtained using the scaling laws (Utsu 2001) from the seismic moment of the F-net moment tensor solution d Values obtained from the residuals in initial fault parameters and the number of data

Fig. 2
Fig. 2 Data fitting for terrestrial geodetic data with the optimal fault model in a the joint inversion, b the GNSS inversion, c the tsunami-source inversion.The black and red arrows indicate the observed and synthetic static horizontal displacements, respectively.The black star denotes the epicenter.The black frame drawn in purple and the thick black line represent the outline of the optimal fault model and the upper edge of its fault, respectively.The blue arrow indicates the slip vector of the optimal fault model.The circles denote the hypocenters of the aftershocks (M ≥ 3) within one day after the mainshock.d-f Enlarged views of the dotted region in a-c.The gray-colored focal mechanism represents the F-net moment tensor solutions of the 2016 off-Fukushima earthquake.The blue-line beachball represents the focal mechanism of the optimal solution for each inversion

Fig. 3
Fig. 3 Spatial distributions of the uplift and subsidence for a Kubota et al. (2021), b the synthetics of the joint inversion, c the synthetics of the GNSS inversion, and d the synthetics of the tsunami-source inversion, and spatial distributions of the residuals of e the joint inversion, f the GNSS inversion, and g the tsunami-source inversion.The star denotes the epicenter

Fig. 4
Fig. 4 Relationship of fault-parameter ensembles in the joint inversion with heat maps of fault parameters and histograms of each fault parameter.The cross in the heatmaps and the dashed line in the histograms indicate the fault parameters of the optimal solution having the maximum posterior probability among solution ensembles.M w , VR GNSS , and VR TS were not included in the model parameters, but calculated from the model parameter values for each ensemble

Fig. 5
Fig. 5 Comparison of fault-parameter relationship among the joint inversion (red, Fig. 4), the GNSS inversion (black, Additional file 1: Figure S1), and the tsunami-source inversion (blue, Additional file 1: Figure S2) shown by the relationship of each parameter pair and histograms of each parameter.Solid contours indicate the outlines of the posterior possibility distributions (areas with a posterior probability of ≥ 0.01%) for each inversion.The cross and dashed line model parameters of the optimal solution for each inversion

Fig. 6
Fig. 6 Spatial distribution of the ensembles of fault models obtained by a the joint inversion, b the GNSS inversion, and c the tsunami-source inversion.The upper panels indicate the spatial distribution of frequency for fault-plane ensembles (color distribution) with the optimal fault model (black frame).The lower panels indicate the selected 40 fault models (gray frame) with the optimal fault model (black frame drawn in purple).The star represents the epicenter.The circles in the upper panels denote the hypocenters of the aftershocks (M ≥ 3) within one day after the mainshock.The gray-colored focal mechanism represents the F-net moment tensor solutions of the 2016 off-Fukushima earthquake.The blue-line beachball represents the focal mechanism of the optimal solution for each inversion

Table 1
Settings of initial values and prior PDFs

Table 2
Optimal fault models in the joint, GNSS, and tsunamisource inversions a Variance reduction of the GNSS data b Variance reduction of the tsunami-source