Numerical experiment toward simultaneous estimation of the spatiotemporal evolution of interseismic fault slips and block motions in southwest Japan

Interseismic strain around southwest Japan can be attributed to slip phenomena on the plate interface and processes inside the continental plate, such as block motions and inland fault slips. Although many crustal deformation analyses have been carried out, the simultaneous estimation of these phenomena and comprehensive discussion regarding the effect on the total strain budget remain topics to be investigated. In this study, we conducted numerical experiment to evaluate the possibility of simultaneously monitoring the spatiotemporal evolution of interseismic fault slips and block motions based on state-space modeling. We aimed to estimate fault slips and block motions using the time series of dense continuous Global Navigation Satellite System sites covering southwest Japan, encompassing 25 years from 1996 onward. We calculated synthetic block motions of the forearc region, monotonically increasing back slips on block boundary faults, and long-term slow slip events in the Tokai region, Kii Channel, and Bungo Channel. Subsequently, we generated the expected synthetic displacement time-series at GNSS stations. Applying Kalman filtering, we successfully estimated the spatiotemporal evolution of block motions, cumulative back slips, multiple slow slip events occurring in different regions. Although the contributions of fault slips and block motions showed slight trade-off, main characteristics of the slip distributions and the direction of block motions were well-recovered. We recovered slow slip events with slips of 5–10 cm or larger. We investigated the estimation uncertainty and separation precision of the unknown parameters using the covariance matrix estimated by Kalman filtering. Focusing on the structure of the non-diagonal component of the covariance matrix, we evaluated the complex effects of site and subfault locations on the estimated slip spatial distribution bias. For instance, the extent of the correlation between subfaults suggested that the three slow slip regions have different tendency of the uncertainties of slip areas, extending toward landside, seaward, or both of them. Our framework enables the comprehensive evaluation of the contributions and uncertainties of various deformation sources covering the entire subduction zone.


Introduction
Interseismic deformation in a plate subduction zone is attributed not only to slip and coupling on a single plate boundary fault, but also to strain partitioning on overriding and subducting plates (e.g., DeMets 1992;McCaffrey 1992;Shen-Tu et al. 1995).Accurately monitoring spatiotemporal distributions of strain partitioning and/or predicting their future evolutions are essential for evaluating the possibility of future large earthquakes.
Southwest Japan is one of the most well-studied regions for strain partitioning.In this region, the northern margin of the Philippine Sea plate subducts along the Nankai Trough, beneath the southeast margin of the Amurian plate at a rate of approximately 6 cm/year (DeMets et al. 2010) (Fig. 1).The main boundary between these two plates has experienced M8 class megathrust earthquakes approximately every 100-200 years (Ando 1975;Ishibashi 2004).In this region, the direction of the plate subduction deviates approximately 30° counterclockwise from the perpendicular direction to the trench axis.
According to geological (e.g., Kanaori 1990;Tsutsumi et al. 1991) and geodetic studies (e.g., Sagiya et al. 2000;Miyazaki and Heki 2001;Tabei et al. 2003), the trench parallel component of this oblique subduction is partly accommodated by a right-lateral motion of the fault system called the Median Tectonic Line and the Niigata-Kobe Tectonic Zone, which divides southwest Japan into the backarc and forearc regions.
Since the mid-1990s, crustal deformation in southwest Japan has been monitored by a continuous Global Navigation Satellite System (GNSS) observation network operated by the Geospatial Information Authority of Japan (Sagiya 2004;Tsuji et al. 2017), called GNSS Earth Observation Network system (GEONET).In addition, over 20 seafloor GNSS-Acoustic (GNSS/A) stations have been installed by the Japan Coast Guard and Nagoya University since the early 2000s (e.g., Tadokoro et al. 2012;Yokota et al. 2018).Using these rich geodetic observations, many estimations of fault slips on the plate interface have been carried out.For example, Sherrill and Johnson (2021) jointly estimated the coseismic slip distribution and after slip history of the 1944 M8.1 Tonankai and 1946 M8.3 Nankai earthquakes using GNSS, leveling, triangulation, and tide gauge observations ranging more than 70 years.Yokota et al. (2016) first combined onshore GNSS and offshore GNSS/A observations to estimate spatial distribution of interplate coupling.Many other studies (e.g., Hirose et al. 1999;Suito and Ozawa 2009;Ozawa et al. 2016;Takagi et al. 2016;Kobayashi 2014Kobayashi , 2017) ) have detected recurrent slow slip events (SSEs) with durations of a few years in the Tokai region (G in Fig. 1), Kii Channel and Bungo Channel (E and F in Fig. 1).SSEs that only last for several days, termed as short-term SSEs have been also reported (e.g., Nishimura et al. 2013;Okada et al. 2022).
Strain partitioning inside the plates has also been investigated.It is often modeled by the rigid motion of small tectonic blocks inside the plate and creep or coupling on their boundary faults (e.g., Matsu'ura et al. 1986;McCaffrey 2002;Meade and Loveless 2009).In addition, internal deformation of the block itself has also been proposed to describe the contribution of unknown structures inside the block (McCaffrey 2005).The combination of these effects to describe interseismic deformation is generally referred to as a block model.Since its first application by Hashimoto et al. (2000), block modeling have been carried out for southwest Japan or the whole of Japan Islands by some previous studies (Wallace et al. 2009;Loveless and Meade 2010;Nishimura et al. 2018;Kimura et al. 2019).The above studies suggested that approximately one-third of the relative plate motion may be attributed to several block boundaries inside the Amurian plate.
The described studies have made great efforts to reveal the spatiotemporal distributions of fault slips and block motions around southwest Japan, with wide ranges in both spatial and temporal scales.However, combining their results to discuss the interaction between these phenomena and their effects on the total slip budget remains an important topic.One of the reasons is that the previous block modeling studies did not analyze displacement time series but instead used an average displacement velocity field spanning more than a few years as the observation.Therefore, temporal variations within the analysis periods were disregarded, and the estimated coupling distribution may have included the contributions of transient slip phenomena.Other reason is that Fig. 1 Sitemap of our analysis.Red squares denote the GEONET sites used for our main analysis.Thick blue lines represent the assumed block boundaries while thin blue lines represent the other block boundaries, suggested by Nishimura et al. (2018).Grey rectangles represent the assumed subfaults.The contour of the plate interface, with an interval of 10 km, is represented by grey dashed lines (Hirose et al. 2008a, b).The red arrow represents the direction of plate subduction given by the Mid-Ocean Ridge Velocity (MORVEL) plate velocity model (DeMets et al. 2010).The geographical names used in this paper are represented on the right side.We also show the known active fault systems that correspond to the block boundaries.The outline of the Chubu block corresponds to the (a) Itoigawa-Shizuoka Tectonic Line, (b) Atotsugawa fault, (c) Biwako-Seigan and Hanaore fault, (d) Arima-Takatsuki fault, and (e) Rokko-Awaji fault.That of the Shikoku block corresponds to the (f ) Median Tectonic Line, (g) Futagawa-Hinagu fault, and (h) Southern Kyushu Shear Zone most fault slip estimations individually focused on events occurring in different areas and periods, using different geodetic data sets and inversion methods.This resulted in various estimation uncertainties for each of the estimated slip distributions.Although Sherrill and Johnson (2021) evaluated the slip budget around the Shikoku region by combining their results with other estimations of SSEs and interplate coupling, it is difficult to show the uncertainty of their interpretation clearly.
Some studies have already highlighted the importance of simultaneous estimation covering the spatiotemporal evolution of different types of slip phenomena.Ochi and Kato (2013) proposed a method to simultaneously estimate interplate coupling and SSEs using displacement time series without excluding interseismic trends.Applying its method for a long-term SSE in the Tokai region (2000)(2001)(2002)(2003)(2004)(2005), they showed that the seismic moment released by the SSE was clearly smaller than that estimated by previous studies that excluded an interseismic trend in advance and defined the SSE as the remained deviation (Suito and Ozawa 2009).They proposed that analyzing SSEs and interseismic rates separately can overlook the contribution of interplate coupling during SSE occurrence, leading to the misinterpretation of their contribution to the total slip budget.
As explained above, a continuous and simultaneous estimation of the spatiotemporal evolution of various types of fault slips and block motions has rarely been attempted, despite the accumulation of dense long-term geodetic observations.Moreover, the real-time monitoring of interseismic fault slip phenomena has also not been achieved.Considering this background, we established the final goal of our research as the construction of a framework for continuous monitoring of the spatiotemporal evolution of fault slips and block motions covering the entire southwest Japan.Toward this final goal, this paper numerically examines the possibility of simultaneously estimating fault slips and block motions based on

Method
Figure 2 illustrates the entire flow of our numerical experiment.We modified the block modeling formulation by the previous studies to apply it to time variable estimations based on Kalman filtering.Then, we prepared certain given values of fault slips and rigid block motions, hereafter we call "true values".To do so, we once carried out estimation of fault slips and block motions based on the real GNSS time series of southwest Japan.Modifying its results, we prepared the true values and calculated the corresponding synthetic displacement time series at the GNSS sites.Finally, we input the synthetic displacement time series into Kalman filtering and examined the recovery of the model parameters.This section provides details on the formulation of the block modeling and the configuration of the block boundary faults.

Observation equation and estimation algorithm by Kalman filtering
We modified the block modeling formulated by Meade and Loveless (2009) and Savage et al. (2001) to use cumulative displacement time series as observations, instead of the displacement velocity field.We assumed the rigid block motions and elastic strain due to fault slips on block boundary faults.We did not include the internal deformation of the block itself because it showed small contribution for the explanation of the observations.First, rigid block motion was approximated by the rotation of a crust along a spherical surface.We defined a rotation rate vector of the j th block ˙ j that describes the direction of the rotation axis and the rotation rate with its direction and length, respectively.The displacement velocity vector at the i th observation site v B,i was given by an outer prod- uct of the rotation rate vector ˙ j = ˙ xj , ˙ yj , ˙ zj T and the i th site position vector r i = x i , y i , z i T as follows: (1) state-space modeling.We then discuss the key topics that need to be addressed to precisely analyze real data.In the following sections, we first describe the formulation of block modeling and its application to Kalman filtering.Subsequently, we explain the procedure of our numerical experiments and present the results.Finally, we discuss the separation precision of model parameters using the estimated covariance matrix.
Note that ˙ j and r i are the vectors under a global Carte- sian XYZ coordinate, with an origin at the center of the Earth.Contrastingly, v B,i is under local EW, NS, and UD (ENU) coordinates, along with a tangent plane at the observation site.The operator P V ,i transforms the XYZ to ENU coordinates, where φ i and θ i are the latitude and longitude of the i th observation site.G B,i is an operator used to calcu- late an outer product of ˙ j and r i .Although it is not repre- sented in the equation, we divided G B,i by the Earth radius, , in the actual estimations to maintain numerical stability.Then, the product of the Earth radius and rotation rate vectors, |r i | • j , will be estimated as unknown whose unit is m/yr.In the following explanation, we display not |r i | • j but j as the obtained results with the unit of /yr.Although Eq. ( 1) connects a single block and a single site on that block, a real analysis involves more than two blocks and sites.Therefore, the operators P V ,i and G B,i for each site were assigned to compose a design matrix that connects the rotation rate vector of each block ˙ j with the displacement velocity at each site v B,i as follows: (2) Here, the subscripts ns and nb represent the total num- ber of sites and blocks, respectively.The position of each of G B,j in the column direction of G B controls which block the site locates.Here, the 1st to nsI th sites locate on blocks and others locate in the exterior fixed region.
Secondly, the elastic deformation due to slip on the block boundary fault was described.Meade and Loveless (2009) introduced a connection between the interseismic slip deficit and the differential block motion of two blocks on either side of that fault.However, our purpose is the continuous estimation of the temporal evolution of both interseismic slip deficit and transient slip.Therefore, the connection between fault slips and block motions is expected to be further complicated.In this study, we did not introduce any connections for the simplicity of the algorithm.The development of the complete formulation of time-variable block modeling is future work.While Meade and Loveless (2009) used a triangular dislocation model suggested by Meade (2007), we assumed rectangular faults and used the elastic half-space model proposed by Okada (1992) for the calculation of Green's function, G s .We allowed free slip direc- tion so that the unknown parameters are the strike and dip slips of each of the l th subfault, ss l and ds l .Although we adopted no specific constraint, the coupling on the fault is expected to be described by back slips in nearly the opposite direction to relative block motion (Savage 1983).Cumulative displacement for the sites is given as follows: Here, the subscript nf represents the number of sub- faults.Green's function G i,ss,l and G i,ds,l give the ENU dis- placement at the i th site generated by unit strike and dip slips on the l th subfault.
Finally, the interseismic cumulative displacement from the initial epoch t 0 to time t , d(t) was estimated as sum of block motions and fault slips: Note the difference in dimension of �(t) and s(t) ; the for- mer represents velocity, while the latter represents cumulative quantity.It should be also noted that �(t) is not a momentary rate at time t but rather the average rate from t 0 to t .In this context, the elapsed time t − t 0 is multiplied by the P V G B to convert velocity into cumulative amount.In case of fault slips, we directly designate the cumulative slip from t 0 to t , s(t) as an unknown parameter to help intuitively understanding the obtained results.We used cumulative displacement as observations because using momentary displacement velocity may result in the amplification of various non-tectonic noise factors.We adopted a spatial smoothing constraint for the slip on each subfault.L represents the Laplacian operator for the smoothing constraint.We applied Eq. ( 4) to state-space modeling, as shown in Eq. ( 5).Hereafter, subscript k indicates the epoch number. (3) where the state vector x k , observation vector d k , and design matrix H k are defined as follows: The first term in Eq. ( 5) represents the state equation.As we did not assume any control physics of the system, the state vector x k was transferred to the next epoch k + 1 , with only the addition of system noise α k .The sec- ond term in Eq. ( 5) represents the observation equation, as described in Eq. ( 4).t and β k represent the sampling interval of the analysis and observation noise, respectively.The system noise covariance matrix Q k and obser- vation noise covariance matrix R k are defined as follows.
where σ B and σ S represent the process noise of the rota- tion rate vectors and fault slips, respectively.We adopted a common value of σ B for all components of rotation rate vectors for all blocks, and a common value of σ S for both strike and dip slips on all subfaults.In this study, (5) , we assumed a white noise stochastic process for all unknowns to handle various time scales and dynamic ranges.σ e , σ n , and σ u represent the observation noise of the EW, NS, and UD components of the displacement time series, respectively.We set their amplitudes as 0.5 cm, 0.5 cm, and 2.5 cm respectively, and adopted homogeneously for all used sites.These are nearly twice as large as the average standard deviation of real data from the GNSS sites around southwest Japan, which is the setting to avoid overfitting of the data.Finally, κ con- trols the strength of spatial smoothing for fault slips.

Settings of block boundary faults
We followed the configuration of block boundary faults by Nishimura et al. (2018).For the inland faults, it was originally based on a combination of known surface fault traces, the spatial distribution of shallow seismicity, and the strain rate observed geodetically by the previous studies.For the subducting plate boundary, the models of Baba et al. (2002) and Hirose et al. (2008a, b) were used.
We assumed two blocks for the east and west halves of the forearc region of southwest Japan, which we hereafter refer to as the "Chubu block" and "Shikoku block, " respectively (Fig. 1).As indicated by (a) to (h) in Fig. 1, the assumed block boundaries basically correspond to the known active fault systems.Strictly, we should also set the block motions of the exterior regions.However, it requires additional observation sites and considerations of block boundary faults, since we independently estimate fault slips and block motions.Then, the analysis region can become larger and larger to cover the entire Japan.The feasibility of applying Kalman filtering in such a large scale of analyzed region, the number of observations, and length of state vector remains uncertain.Therefore, we treated the exterior regions as fixed.This simplification should be discussed in the future.We set 105 rectangular subfaults surrounding these two blocks (Fig. 1).For all sections, we placed three subfaults in the dip direction, with a depth range of 0.5 to 24.5 km and a dip angle of 30° to nearly vertical.The typical length and width of the subfaults are 10-50 km and 8-16 km, respectively.Hereafter, these are referred to as "inland subfaults, " and free slip direction were allowed for them to avoid kinematic inconsistency between adjacent fault segments with different strikes.We assumed creeping for the boundary between the Chubu and Shikoku blocks.For the plate interface, we set 578 subfaults from the trench axis of the Nankai Trough to a depth of 60 km, with length and depth of mostly 20-50 km and 5-70 km.These subfaults are hereafter referred to as "interplate subfaults, " for which we also assumed free slip direction to describe strain partitioning from the direction of plate subduction.For both inland and interplate subfaults, their lengths were set to be sufficiently larger than the average interval of the GEONET sites.In this study, as an antonym of "back slip", we call reverse slips on the interplate subfaults as "forward slip".In total, we considered the cumulative time series of strike and dip slips for each 683 subfault.Adding the rotation vectors of the two blocks, the length of the state vector to be estimated was 1372.

Numerical experiment
In this section, we discuss the preparation of "true values" of fault slips and block motions for the input of the numerical experiment.It is desirable that these true values precisely approximate the real tectonic process occurring in southwest Japan.However, the simultaneous estimations of fault slips and block motions across the entire southwest Japan have been rarely conducted except for Nishimura et al. (2018).Projection of their results into our analysis is not easy and can cause various arbitrariness due to significant differences in block boundary settings and subfault configurations.Moreover, our analysis requires a time series of model parameters, not just a snapshot at a specific time.On the other hand, one of the largest problem in time-variable estimation in Japan is the postseismic effect of the 2011 Tohoku-Oki earthquake.To check its impact and discuss future topics toward real data analysis, examining the response of Kalman filtering toward real data is beneficial.
Based on the above motivations, we conducted an estimation of fault slips and block motions using real GNSS time series.We then modified its results to prepare the true values and generated corresponding synthetic displacement time series.This stage also involved the selection of hyperparameters.

Pre-estimations of fault slips and block motions
We simultaneously estimated the temporal evolution of fault slips and rigid block motions from daily coordinate time series at the 743 GEONET sites (Fig. 1).We obtained time series from 1996 to 2021 using precise point positioning with ambiguity resolution analysis on the Gip-syX-1.4software.We set the site 0462 (Fukue) as a fixed site and removed the offsets due to site maintenances and coseismic displacements resulting from known large earthquakes.We did not remove any postseismic effects of large earthquakes and seasonal fluctuation.
We set a process noise value of σ S = 0.5 m and adopted a white noise stochastic process to handle slip phenomena with various velocities.We set the spatial smoothing parameter κ as 0.04 for interplate subfaults and 0.02 for inland subfaults, respectively, based on the trade-off curve between displacement fitness and smoothness of the obtained slip distribution (see also Additional file 1: Figs.S1 and S2 for details).Since the inland subfaults are much closer to the surface, they are sensitive to the motion of a few neighboring observation sites.Therefore, to suppress local overfitting, we adopted stronger smoothing for the inland subfaults.For the rotation rate vectors, we set the process noise value of σ B = 1.57e−9/ yr (0.01 m/yr), based on the typical order of magnitude of the rotation rates estimated by Nishimura et al. (2018).We consider that the temporal fluctuation of rotation rate is very small or difficult to detect in reality.However, to catch up the burn-in behavior in the earlier phase of the estimation, we also adopted white noise stochastic process.The above hyperparameter settings will be also applied to the main estimation.
Figure 3a shows the estimated slip distribution on March 4, 2011, 1 week before the Tohoku-Oki earthquake.The locations of strong couplings are well consistent with previous studies (Nishimura et al. 2018;Yokota et al. 2016), that were estimated around the Tokai to Tonankai region, nearby Kii Channel, Nankai region, and Bungo Channel.We also observe some unrealistic large slips in the downdip subfaults beneath the Chubu, Kinki, and southern Kyushu regions (A and B in Fig. 3a).They reproduce the westward movement of backarc regions outside of the assumed fault region and the rotation of the Kyushu region, respectively.Considering these characteristics, we only adopted the slips on the inner part of the assumed fault region for the "true slips", as shown in the lower panel of Fig. 3a.We excluded slips with depths shallower than 10 km and deeper than 40 km, corresponding to the shallowest one row and deepest five rows of subfaults.Similarly, we excluded the subfaults of the eastern three columns and the western four columns.The inland subfaults generally showed a left-lateral slip which describes the coupling against the right-lateral movement of the forearc regions.Subfaults in the Itoigawa-Shizuoka Tectonic Line showed unstable large slips owing to imposed westward movements due to the subduction of Pacific plate (C in Fig. 3a).Therefore, we excluded all slips of this section.Other sections partly show large depth variation of the slips due to poor resolution along dip caused by higher dip angle.Then, we finally averaged the slips on each of the three subfaults in the dip direction to remove the depth variations.The resulting spatial distribution of the assumed true slips is very similar to that estimated by Nishimura et al. (2018).The amount of the true inland slips is generally 2-17 mm/yr.Strictly speaking, the directions of the slip deficit and the differential block motion should have consistency as we mentioned in the previous chapter.In case of our pre-estimation results, most slips showed the sense of back slip which is opposite to the differential block motions.However, the block motions generally direct more west or southwest, carrying nearly pure lateral motions which is parallel to the trench.The difference of azimuth reached 40°-80° in the interplate subfaults around the Kii Channel, exceeding the estimated standard deviation of fault slips.In case of the inland subfaults, significant dip slips were estimated whereas the differential block motion mostly corresponds to strike slip.Further consideration will be required for the physical adequacy of the "true values" and the estimation algorithm.
The estimated slip time series generally shows a monotonically increase except for the known long-term SSEs (Fig. 3b).In addition, slip time series show a log-like curve reflecting the postseismic effect of the 2011 Tohoku-Oki earthquake.Deviations of the slip time series from the pre-seismic trend reached approximately 0.1-0.2m for the main coupled region.Figure 3c shows the estimated time series of the rotation rate vectors.They all showed larger fluctuations in the first 5 to 10 years, after which they gradually converged since the signal-to-noise ratio of the estimations gradually increases.We can also observe trend change at the time of the 2011 Tohoku-Oki earthquake and 2016 Kumamoto earthquake, reflecting its postseismic effects.The estimated rotation rate vectors on March 4 give displacement velocity of approximately 8-16 mm/yr for the Chubu block and 3-4 mm/ yr for the Shikoku block respectively, that are similar or smaller than Nishimura et al. (2018).The displacement by rigid block motions mainly show the direction from WNW to WSW (Fig. 4b), describing the lateral motion of the forearc region which is usually referred to as sliver motion (e.g., McCaffrey 1992).In total, our pre-estimations reproduced the observations with a root mean square (RMS) of 6-7 cm (Fig. 4d).

Generation of the "true values" and synthetic displacement time series
As described above, slip time series estimated from the real data generally show a monotonically increase and large postseismic effect of the Tohoku-Oki earthquake.Rotation rate vectors showed burn-in behavior due to the gradual improvement of the signal-to-noise ratio and then converged to certain values.According to such results, we did not directly use the entire pre-estimated time series but instead used the values on March 4, 2011 as the initial reference values, ˙ ref0 and s ref0 .Using them, the time series of the true values were given by Eq. ( 8).
We assumed constant values for the true rotation rate vectors throughout the analysis period.For the true ( 8)  S3 for details).We added a flat ramp whilst maintaining the final slip values given by Eq. ( 8).We defined the slip of SSEs as the loss of cumulative back slip resulting from SSE occurrence, hereafter referred to as "relative SSE slip." As this depends on the long-term back slip velocity for each subfault, we also checked the apparent slip variation during the SSE period in the numerical experiments, hereafter referred to as "apparent SSE slip." We set the slip regions, slip amounts, and periods of the SSEs roughly based on previous studies (e.g., Ohta et al. 2004;Miyazaki et al. 2006;Liu et al. 2010;Ozawa et al. 2013;Kobayashi 2017), although they can contain slip overestimations as we mentioned in the introduction part.For the Tokai region, we assumed two events with durations of 4.5 years with a common slip region and a maximum slip of 30 cm for both events.Similarly, for the Bungo Channel, we assumed three events with durations of 1.5 years and a maximum slip of 27 cm.The corresponding moment release of the five events was around a moment magnitude (M w ) of 7.0.They give an approximately 10 cm of surface displacement in maximum.Contrastingly, we assumed much smaller slip amounts for SSEs in the Kii Channel.We assumed two events lasting 2.5 years and 3 years, with maximum slips of 2 cm and 5 cm respectively.The corresponding moment release was M w 6.4-6.6.They give 1.0-1.5 cm of maximum surface displacement, which is only slightly larger than the noise level of the displacement time series.
We generated the synthetic displacement time series from the time series of true values using Eq. ( 9).We simply added displacement by true fault slips and rotation rate vectors and set no specific weight between them.We added the synthetic observation noise β syn (t) with a power inversely proportional to its frequency, which is referred to as pink noise or flicker noise.We set noise amplitudes of 0.26, 0.25, and 1.25 cm for the EW, NS, and UD components, respectively.These are the average standard deviations of the real displacement time series at the 743 sites used in the pre-estimation. (9)

Results of numerical experiment
Using the generated synthetic displacement time series as input observations, we carried out the simultaneous estimation of the fault slips and rotation rate vectors to examine the recovery of them.Figure 5 compares the estimated and true back slip distributions on March 4, 2011, approximately 14.5 years from the beginning of analysis (see also Fig. 1 for the geographical names used in the explanations).The cumulative back slip is smoothed for the shallow subfaults due to the spatial smoothing constraint and poor resolution of the GNSS for the offshore regions.In contrast, slips on downdip subfaults at depths greater than 15 km, including those in the region of long-term SSE occurrence, were generally well-recovered within 0.1-0.2m and 10-20% of residual.Peaks of back slip around Nankai and Bungo Channel show the best recovery.Although the above slip residuals still exceed the standard deviation of fault slips estimated by Kalman filtering, they are small considering that they are cumulated amounts for more than 14 years.Back slips on the inland subfaults mostly show residual slips within 2-5 cm and comparable to the estimated standard deviations.Thus, we conclude that back slips on the inland subfaults can be well-recovered.
Table 1 shows the estimated and true values of the rotation rate vectors on March 4, 2011.Although the norm of the vectors showed slight overestimation or underestimation, their directions well-agreed.The estimated time series required 4-5 years to converge into the true values.Thus, when real data is analyzed, rotation rate vectors should be set as temporally changing variables to accommodate for such burn-in behavior, even in cases that they are truly constant.
We show residuals between the estimated and true displacement fields in Fig. 6.The total residual displacements show no clear tendency for their directions (Fig. 6c).RMS for the sites in the forearc two blocks are around 6 mm for the horizontal component.However, the residual displacements resulting from fault slips were mainly oriented west to southwest (Fig. 6a), with RMS of 10-11 mm.Oppositely, those due to rigid block motions were oriented east to northeast (Fig. 6b), with RMS of 5-6 mm.These opposite directions indicate that a trade-off between fault slips and rigid block motions occurred, associated with the sliver motion of forearc region.
Next, we focus on the estimations of the SSE slips.Figure 7 compares the estimated and true relative slips of the SSEs (see Additional file 1: Fig. S4 for a comparison of the apparent slips).In general, the SSEs in the Tokai region and Bungo Channel were well-recovered.The relative slips and equivalent released moment showed underestimations of approximately 15-20% and 10-20% respectively, partly corresponding to the underestimations of back slips.The slip recovery are similar for second and third events.In the case of the SSEs in the Kii Channel, it seems much more difficult to recover the fault slips within just a few cm.Slip distributions seem to be unstable and are clearly extended toward the updip side of the true slip area.Such results are reasonable because the estimated standard deviation of the fault slip around the Kii Channel is approximately 3-4 cm, as shown in Fig. 10a in the subsequent discussion part.Thus, the 2 cm slip in the first event is completely below the uncertainty and the 5 cm slip in the second event is slightly above the uncertainty.From these results, we estimate that the detection limit of our analysis is an SSE with slips of at least 5 cm or more in 1 to 5 years.Then, we can precisely estimate the slip distribution in case of more than 10 cm slips.
Examples of the estimated and true slip time series are shown in Fig. 8.The ratio between the estimated and true slip values was generally constant, except for the first 5 years, indicating that the estimation values deviated from the true values in a non-gradual yet constant manner throughout the analysis.Notably, we observed some leakage of SSE slips toward the outside of the true slip areas.For the Tokai region, slight transient changes mainly appeared on the downdip side of the true slip area.The ratio between the true and estimated slips shows slight bends that correspond to the SSE occurrence periods, although the step changes of the estimated slip time series itself is not so clear (dashed black lines in Fig. 8b).In case of the Bungo Channel, the clearest slip leakages occur on the west of the true slip area (Fig. 8d).In addition, both the updip and downdip side also show some slip leakages.Finally, the Kii Channel shows slip leakage mainly toward the updip side (Fig. 8f ).As described above, three SSE regions show different tendency of estimation uncertainty of the slip area, extending toward the seaward, landside, or both of them.In summary, we were able to successfully estimate the spatiotemporal evolution of cumulative back slips, multiple SSEs occurring in different regions, and rigid block motions covering the entire subduction zone.The spatial distribution of the cumulative back slips was well-recovered except for the offshore region nearby trench.In addition, the directions of the rigid block motions were precisely recovered.However, we observed a slight tradeoff between fault slips and rigid block motions associated with the explanation of the forearc sliver motion.For transient slips, we were able to detect long-term SSEs with slips of at least 5 cm or more.We observed slip leakage toward the outside of the true slip regions.Then, the three regions showed different tendencies of leakage.

Separation precision of the unknown parameters evaluated with a covariance matrix
Even though there are many fault slip analyses, few studies have discussed the structure of their covariance matrix (e.g., Bletery et al. 2016).This section discusses the estimation uncertainty and separation precision of the unknown parameters through the structure of a covariance matrix estimated inside Kalman filtering.Figure 9 shows the correlation coefficient matrix that was converted from the covariance matrix on the end of the analysis period.We observed almost no temporal change in their values, indicating that the covariance matrix is mostly determined by the initial model settings and does not depend on the characteristics of the input observations.The indices of the unknown parameters are in the following order: interplate subfaults, inland subfaults, and rotation rate vectors.The subfaults are ordered from shallower and eastern to deeper and western.Therefore, the diagonal stripes in Fig. 9 represent the correlation between subfaults in the dip direction, while the thickness of each stripe describes the correlation in the strike direction.The gradual decrease in the correlation from the center to the bottom right side of the matrix indicates a weaker correlation between distant subfaults, which is mainly reflecting the effect of spatial smoothing.The vertical pattern present in the right-side end shows the correlation between fault slips and rotation rate vectors.
To clearly visualize the overall characteristics of slip uncertainty, we generated additional two types of figures (Fig. 10).Although we only show the results for strike slips, those for dip slips are very similar (see Additional file 1: Fig. S5). Figure 10a shows the estimated standard deviation of the fault slips extracted Its spatial pattern is similar to that of the diagonal component of the resolution matrix calculated from the Green's function (Fig. 10b).The correlation coefficient between these two maps is around 0.8, indicating that the standard deviation mainly reflects the resolution of slip existence or its amount, that are controlled by the responsibility toward surface displacements of the nearest observation sites.Indeed, standard deviations are larger in the offshore and downdip regions that are distant from GNSS sites.Subfaults with higher dip angle especially show larger standard deviation, since they are difficult to generate large horizontal displacement by dip slips.
Figure 10c shows the total amount of correlation between each subfault and all other subfaults, hereafter referred to as the "correlation area." The correlation area of the i th subfault was calculated using Eq. ( 10). (10 Here, CC ij represents the correlation coefficient between the i th and j th subfaults, obtained from the non-diagonal component of the covariance matrix.S j and nf represent the area of the j th subfault and the total number of subfaults, respectively.Therefore, Eq. ( 10) provides the sum of the products of the absolute correlation coefficient and area of each subfault, excluding the selected subfault itself.Correlation area may describe the uncertainty of the area of slip centered at each subfault.It is not a simple resolution but a more complex effect caused by locations and layouts of subfaults, blocks, and sites.In addition, the impact of the spatial smoothing constraint is also included.We discuss the interpretation of correlation area in the next section.
The spatial patterns shown in Fig. 10a, c clearly differed with correlation coefficient below 0.5-0.7.For instance, the standard deviations in the Tokai region and Kii Channel show similar amplitudes, whereas those in the Bungo Channel are approximately 1.2 times larger.In case of the correlation area, the Tokai region shows the smallest values and other two regions show 1.5-2.0times larger values.These characteristics suggest that the diagonal and non-diagonal components of the covariance matrix provide different indexes of slip uncertainty.We were able to comprehensively evaluate the slip uncertainty of different regions in the entire subduction zone.This is the main advantage of our simultaneous estimation covering the entire subduction zone.
In Fig. 11, we show the extent of the correlation between the three SSE areas and the surrounding subfaults.For each of the surrounding subfaults, we considered multiple correlation coefficients associated with multiple subfaults within the focused area.Therefore, we took their average absolute value as a representative.In the Tokai region, relatively strong correlations were observed mainly on the downdip sides of the SSE area.This probably indicates that the estimated slip areas tend to be drawn landward, which is closer to the GNSS sites and has higher sensitivity.Contrastingly, the slips in the Bungo Channel show strong correlations mainly on the west and updip sides, indicating that the estimated SSE slips tend to extend seaward with lower resolution.The downdip side also shows relatively strong correlations.Finally, the Kii Channel shows a strong correlation on the updip and downdip sides.These tendencies are consistent with our numerical experiment results.They suggest slip leakage toward the downdip side for the Tokai region (Fig. 8b), toward the west, updip, and downdip sides for the Bungo Channel (Fig. 8d), and toward the updip side for the Kii Channel (Fig. 8f ).As shown above, the three SSE regions show different tendencies of slip area uncertainties.It can then be evaluated by the structure of the covariance matrix.

Interpretation of "the uncertainty of slip area" and its future use
In the above discussion, "the uncertainty of the slip area" includes the combined effect of the physics of the crust described by Green's function and the layouts of the sites and block boundaries.It is then further filtered by the effect of the spatial smoothing constraint and the subfault setting.For instance, larger subfaults and stronger spatial smoothing make the correlation areas larger.The above effects, which constitute "slip area uncertainty" should be precisely separated in the future.
In such a situation, it is difficult to determine a specific absolute value of the uncertainty of the slip areas.On the other hand, we have confirmed that the relative characteristics of the correlation areas and correlation extents are generally stable under different smoothing strengths (see Additional file 1: Fig. S6).While correlation areas are overestimated under too strong smoothing, their values reach certain minimum under sufficiently weaker smoothing strengths.Presumably, these are the real values of the correlation areas, which are controlled by the Green's functions and the model geometry.Therefore, one idea is to adopt the smoothing strength at the corner point of such an L-shaped curve as the optimum.However, this idea assumes that the Laplacian smoothing is the best form to describe the slip distribution.Some previous studies (e.g., Mai and Beroza 2002) have statistically analysed a number of slip estimations and reported the existence of spatial correlations in real fault slip.Consideration of the smoothing method is an important future topic.
The method of subfault discretization is closely related to the smoothing constraint.We expect that the real fault slips do not have sharp discontinuities, like the "true SSE slips" we have assumed, which cause infinite shear stress.On the other hand, it is difficult to distinguish slip structures at such a small scale and discuss their physical adequacy based on the GEONET sites with an average interval of only 20-30 km.Therefore, a balance between the information contained in the observations and the complexity of the model is essential (e.g., Xu 2022).Although it is beyond the scope of this study, the mathematical determination of the optimal number and size of subfaults is an important issue.For example, Agata et al. (2022) selected the size of subfaults based on the widely applicable Bayesian information criterion (Watanabe 2013).
Another important issue is the introduction of model uncertainties, such as the inaccuracy of Green's function and block boundary settings.Murakami et al. ( 2022) calculated Green's functions based on threedimensional crustal structure model for southwest Japan with systematic changes in the physical properties (density and velocity) of the crust and mantle, to evaluate their effects on slip inversion.They concluded that changing physical properties gives no significant effects exceeding the observation error levels in case of interseismic slip deficit inversion, probably because of the lower S/N ratio and the effect of regularization.Based on these results, we have omitted the introduction of model error effects in the system covariance matrix.However, unlike their experiments, we use cumulative displacements for more than 20 years, which means that the observed amplitude is larger and some effects may still appear.There have been many efforts to construct accurate Green's functions for the Japanese Islands, including viscoelastic effects (e.g., Ichimura et al. 2016;Hori et al. 2021).Finally, we should introduce such results and try to adopt various boundary settings to evaluate their effects.Since the plate interface beneath southwest Japan has sharp bends, rectangular subfaults introduce singularities in the regions of gaps and overlaps.In this study, the subfaults around the Kii Channel and the downdip region below the Kinki region have the largest overlap, reaching 30% or more, which may cause attenuation of the estimated slip amounts.It is better to check the effect of such overlaps or to introduce a more contiguous geometry like triangular mesh (e.g., Meade 2007).
Another important simplification in our analysis is that we have assumed constant stochastic process and process noise values for all unknowns in our analysis.However, time-variable stochastic processes such as Fukuda et al. (2004) can be also adopted to analyse slip phenomena with various time scales and amplitudes.
As mentioned above, the main purpose of this paper is not to construct fully sophisticated models and realistic slip inversions, but to consider a methodology for simultaneous analysis of the entire subduction zone.Our results are not the ground truth, but just the results for the current settings.In the future, the evaluation of slip uncertainties shown in this paper should be repeated with further refinement of the model construction.Using the correlation area and/or extent for different slip areas, as shown in Figs. 10 and 11, an ensemble of the expected bias of different slip areas is generated.By comparing this with the slip distributions estimated from the real data, we can assess whether or not the slip on each subfault is true or not.This can also be applied to evaluate the impact of the joint use of different geodetic observations and/or the selection of the location of newly installed sites.For example, we can search for site locations that minimize the correlation area.Kalman filtering originally allows the combination of multiple observations with different time intervals.Therefore, it would be possible to combine GNSS with other geodetic observations, such as InSAR (e.g., Bekaert et al. 2016;Aslan et al. 2019) and GNSS/A (e.g., Yokota et al. 2016).Time-variable joint inversion has not been performed as much and is also an important future topic.

Other topics toward real data analysis
The synthetic displacement time series used in our numerical experiments do not contain the entire complexity of the real data.For instance, as we have checked through the pre-estimation, the postseismic effect of the 2011 Tohoku-Oki earthquake is the biggest problem of the real GNSS time series around the Japan.If we analyse the real data for more than 25 years including after 2011, the postseismic effect must be removed in some way, such as using the previously suggested deformation models (e.g., Sun and Wang 2015).
In this study, we used the cumulative displacement since June 1, 1996 as an observation.However, more than half of the used sites were established after this reference day.Therefore, we had to estimate the displacement from June 1, 1996 to the day of site establishment and add it to the displacement time series.We attempted spatial interpolation of the displacements of neighboring older sites and temporal extrapolation of the displacement time series of the focused site itself.With these two approaches, we concluded that the initial offset had at least a few centimeters of uncertainty and provided at least a few percent of bias to the estimated cumulative back slips.Such treatment of the heterogeneity of observations is an important point for long-term analysis.
Although we added pink noise with a common amplitude to the synthetic displacement, the noise characteristics of the real GNSS time series are further complex.The power spectrum of real GNSS time series does not seem monotonous but instead bends at longer periods (e.g., Ohta 2016).Thus, the gradient of the spectrum differ from site to site.We can further evaluate the effect of the spatiotemporal noise variation on the estimation results or incorporate it into the structure of the observation covariance matrix.Alternatively, we can also add noise factors such as seasonal fluctuations into the unknown to estimate simultaneously with true deformation sources.
With addressing the above topics, we expect our framework will lead to continuous fault slip monitoring to predict future large earthquakes.In case of real-time monitoring, it can be also used for forecasting subsequent earthquakes just after the first large earthquake, which is especially concerned in the future megathrust earthquake in southwest Japan.

Conclusions
In this study, we considered a method to continuously monitor the spatiotemporal evolution of interseismic fault slips and block motions covering the entire subduction zone in southwest Japan based on state-space modeling.We carried out numerical experiment in the southwest Japan subduction zone, to examine the possibility of using such a framework.We modified the formulation of block modeling developed in previous studies to apply Kalman filtering.We assumed block motions of the forearc region of southwest Japan together with fault slips on 683 small subfaults covering the plate interface and inland faults surrounding the blocks.We attempted the estimation of these model parameters using the time series of 743 continuous GNSS sites covering southwest Japan for over 25 years.As the true fault slips, we set monotonically increasing back slips and seven long-term SSEs occurring in the Tokai region, Kii Channel, and Bungo Channel.For rigid block motions, we set motions toward the west corresponding to the forearc sliver motion.
From these true values, we calculated the expected cumulative displacement time series of 743 GNSS sites as synthetic displacement time series.As a result of their input into Kalman filtering, we confirmed that the temporal evolution of cumulative back slips, SSEs occurring in the three regions, and block motions could be successfully estimated.The total contributions of fault slips and block motions show slight trade-off in the explanation of forearc sliver motion.However, main characteristics of the slip distributions and the direction of rigid block motions were well-recovered.We could estimate the long-term SSEs with slips of 5-10 cm or more.Using the covariance matrix estimated by Kalman filtering, we evaluated the estimation uncertainty and separation precision of the unknowns.By focusing on the structure of the non-diagonal component, we were able to evaluate the complex effects of subfault locations and layouts.For example, the extent of the correlation between subfaults suggested that the three SSE regions have different tendencies of the uncertainties of slip areas.Slips in the Tokai region, Bungo Channel, and Kii Channel tend to extend toward the landside, seaward, and both of them respectively.Our framework enables a comprehensive evaluation of the contributions of various deformation sources among long-term and wide-area interseismic strain.Moreover, our numerical experiment provide helpful information regarding the uncertainties of slip distributions, enhancing the interpretation of real data analysis.

Fig. 2
Fig. 2 Flow chart of the numerical experiment

Fig. 3
Fig. 3 Results of the pre-estimations.a Cumulative slip distribution on March 4, 2011 and assumed "true slip distribution".Slip amount is denoted by a color scale with different color scales for the interplate and inland subfaults.Areas A-C, encircled with dashed lines, indicate unrealistic slips which we excluded when preparing the true slips.b Examples of the estimated cumulative slip time series for the subfaults in the Tokai region and Bungo Channel.The increase of back slip was set in the positive direction, and three subfaults were chosen along the strike direction, denoted by colored rectangles in panel a.Colors of the rectangles and lines are corresponding.Dashed black lines indicate transient variations corresponding to the known long-term SSEs that occurred in each region.c Estimated time series of the rotation rate vectors.Light red, blue, green, and grey lines denote the results of the Chubu block, and dark color lines denote those of the Shikoku block.Black and grey lines show the total rotation rate calculated from the three components of the vectors

Fig. 4
Fig. 4 Horizontal displacement fields calculated from the pre-estimation results.Displacements a by fault slips, b by rigid block motions, c their sum, and d residual between the observations and calculations.Scale of arrows differ for a-d.Note that the figures show displacement not from the "true values" adopted for the numerical experiment but from pre-estimated values

Fig. 5
Fig. 5 Comparison between the true and estimated slip distributions on March 4, 2011.a True and b estimated slip distributions.c Residual between the true and estimated slip distributions, which is the estimated minus true.Red colors represent larger slip for all panels.Note that the color scales for a-c are different.For the visibility, we inserted the horizontal view for the Arima-Takatsuki and Rokko-Awaji fault (d and e in Fig. 1), whose dip angle is nearly vertical

Fig. 6
Fig. 6 Residual between the true and estimated displacement fields.We show the results for a displacements by fault slips, b that by rigid block motions, and c total displacements.RMS values for the horizontal component are written in each panel

Fig. 8
Fig. 8 Comparison between the estimated and true slip time series.The figures show six subfault examples marked by grey or black rectangles in Fig. 7. Figures a-f show the subfaults indicated by a-f in Fig. 7, respectively.Black and red lines represent the true and estimated time series, respectively, corresponding to the left-side vertical axis.Light red lines represent the ratio between the true and estimated slip values, corresponding to the right-side vertical axis.Horizontal, grey dashed lines denote the ratio of one, indicating that the estimated and true slip values are the same.The range of the vertical axes is different for each figure.Black dashed rectangles indicate the slip leakages mentioned in the main text

Fig. 9
Fig. 9 Correlation coefficient matrix on August 9, 2021, converted from the estimated covariance matrix.The structure of the matrix and order of the index of unknown parameters are presented on the left panels

Fig. 10
Fig. 10 Slip estimation uncertainties as described by the covariance matrix, for strike slips.See Additional file 1: Fig. S5 for dip slips.a Estimated standard deviation and b correlation area of the fault slip on each subfault.c Value of the diagonal part of the resolution matrix calculated from the Green's function.The color bar for the resolution matrix and correlation area is not linear.Black lines indicate the assumed slip area of long-term SSEs

Fig. 11
Fig. 11 Spatial extent of correlation between three SSE areas and all other subfaults.From the results for strike and dip slips, we show dominant component of the SSEs in each region.Each panel show the results for a dip slips for the Bungo Channel, b strike slips for the Kii Channel, and c strike slips for the Tokai region.The mean absolute values of correlation coefficients between each subfault and the multiple subfaults inside the focused regions were considered

Table 1
Estimated values of the rotation rate vectors(March 4, 2011)