Data assimilation for fault slip monitoring and short-term prediction of spatio-temporal evolution of slow slip events: application to the 2010 long-term slow slip event in the Bungo Channel, Japan

Monitoring and predicting fault slip behaviors in subduction zones is essential for understanding earthquake cycles and assessing future earthquake potential. We developed a data assimilation method for fault slip monitoring and the short-term prediction of slow slip events, and applied to the 2010 Bungo Channel slow slip event in south-west Japan. The observed geodetic data were quantitatively explained using a physics-based model with data assimilation. We investigated short-term predictability by assimilating observation data within limited periods. Without prior constraints on fault slip style, observations solely during slip acceleration predicted the occurrence of a fast slip; however, the inclusion of slip deceleration data successfully predicted a slow transient slip. With prior constraints to exclude unstable slip, the assimilation of data after slow slip event occurrence also predicted a slow transient slip. This study provides a tool using data assimilation for fault slip monitoring and prediction based on real observation data.


Introduction
Monitoring fault slip behaviors in subduction zones and predicting their short-term evolution are critical for understanding seismic cycles.Specifically, we define the term "monitoring" as estimating the spatio-temporal evolution of fault slips based on various observation data and "prediction" as using physics-based numerical simulations to assess how current slip behaviors might evolve in the imminent future, especially focusing on evaluating spatio-temporal time-series evolution on crustal deformation or fault slips the ongoing slip event.One direct approach to address both tasks is to employ data assimilation (DA) techniques (e.g., Lewis et al. 2006;Fletcher 2022).DA combines observation data with numerical models to obtain a more realistic model that quantitatively explains the observations.DA is widely used in meteorology and oceanography, particularly in practical applications, such as weather forecasting.
DA techniques have been developed and adopted for fault slip estimations in plate subduction zones (Fukuda et al. 2009;Kano et al. 2013;Hirahara and Nishikiori 2019;van Dinther et al. 2019).These previous studies have mostly focused on the optimization problem of frictional properties, which determine the fault slip behavior in physics-based models, and the initial values of simulation variables such as slip velocities on the fault.Through optimization, observation data can be quantitatively explained using physics-based models.Various DA methods have been applied to the fault slip problem, including the Markov chain Monte Carlo (MCMC) method (Fukuda et al. 2009), particle filter (Mitsui et al. 2010;Hori et al. 2014;Banerjee et al. 2023;Marsman et al. 2024), ensemble Kalman filter (Hirahara and Nishikiori 2019;van Dinther et al. 2019;Diab-Montero et al. 2023), and adjoint method (Kano et al. 2013).Most of these previous studies have validated their proposed methods through numerical experiments using synthetic observations, with only a few studies utilizing real observations.DA employing real observations was successfully applied to the afterslip that followed the 2003 magnitude 8.0 Tokachi-oki earthquake that occurred off the Tokachi region, northeast Japan.Fukuda et al. (2009) applied the MCMC method to successfully estimate the posterior probability density functions (PDFs) of frictional parameters, using a single spring-slider system with a sub-daily Global Navigation Satellite System (GNSS) time series for 5 h following the mainshock.Kano et al. (2015Kano et al. ( , 2020) ) developed the adjoint DA method to optimize spatially variable frictional parameters and initial slip velocities from daily GNSS data for 15 days following the mainshock.They also examined the short-term predictability of the afterslip for another 15 days by comparing the misfit between observed and numerically predicted crustal deformations.
Slow slip events (SSEs) are another research target of DA (Hirahara and Nishikiori 2019;Diab-Montero et al. 2023).Many observational studies have inferred that SSEs occurred prior to megathrust earthquakes (e.g., Ito et al. 2013;Ruiz et al. 2014;Radiguet et al. 2016;Voss et al. 2018); therefore, monitoring and predicting the spatio-temporal evolution of SSEs is important, particularly for assessing future earthquake potential.Hirahara and Nishikiori (2019) (hereafter referred to as HN19) first introduced the use of DA for the fault slip estimation problem of SSEs focusing on long-term SSEs in the Bungo Channel, southwest Japan.Through numerical experiments, they demonstrated that the ensemble Kalman filter successfully estimated the frictional parameters on the fault where the SSE occurred.Following HN19, Ohtani et al. (2023) modified the adjointbased method, using a fault model similar to that used by HN19 together with the knowledge that SSEs occur repeatedly, to predict the timing of the next SSE.They numerically demonstrated that the proposed method successfully predicted the timing of the next SSE.The timescale of this prediction was several years considering the recurrence cycle of SSEs.
To our knowledge, this study is the first to apply DA to real SSE observations, in contrast to previous studies.Assuming the fault model in HN19, we assimilated GNSS data recorded during the 2010 Bungo Channel SSE using the MCMC method.Although this fault model is relatively simple, we focus on how such a simple model can explain the GNSS observations and predict their shortterm temporal evolution.We also attempted to predict the temporal evolution during a single SSE, assuming that an SSE is currently ongoing.We aimed to evaluate the short-term evolution of fault slip.This corresponds to the short-term time-series prediction of SSEs with timescales ranging from weeks to a few years at most, which is shorter than the timescale for predicting the occurrence of the next SSE presented by Ohtani et al. (2023).Such short-term prediction will provide information on future temporal evolution of stress transfer in the neighboring region, for example, where coseismic slip will potentially occur.
In the following sections, we first describe the DA settings.We then present the monitoring and prediction results for the 2010 Bungo Channel SSE.Finally, we discuss the predictability of the SSE's short-term evolution.Through this analysis, the present study proposes a strategy for real-time fault slip monitoring and short-term prediction using DA.

Physics-based model of the long-term Bungo Channel SSE
The Bungo Channel is located southwest of Japan in the Nankai subduction zone, where the Philippine Sea Plate subducts beneath the Amurian Plate (Fig. 1a).Beneath the Bungo Channel, SSEs with a moment magnitude of  Hirahara and Nishikiori (2019).This study assumes uniform frictional parameters both inside and outside the SSE patch indicated in the red circle.The locations of GNSS stations are marked by green triangles, with the time series for stations A-D presented in b.The temporal slip evolution, denoted by the dotted line, is shown in Fig. 2d. a Bottom: the tectonic setting of the study area, in which the Philippine Sea Plate (PH) subducts beneath the Amurian Plate (AM) along the Nankai Trough.b Examples of GNSS time series at stations identified in a.The displacements in trench-parallel (X), trench-perpendicular (Y), and vertical (Z) components are represented by green, blue, and black dots, respectively.The red lines depict the calculated time series, derived from 100 frictional parameters sampled from the posterior PDF of the model parameters.Note that the red lines largely overlap, making it visually difficult to distinguish the 100 time series in each component.GNSS = Global Navigation Satellite System ~ 7 were repeatedly observed around 1980, 1985, 1991, 1997, 2003, 2010, and 2018 based on geodetic observations (Hirose et al. 1999;Kobayashi and Yamamoto 2011;Ozawa et al. 2013;Yoshioka et al. 2015;Seshimo and Yoshioka 2022).HN19 introduced a fault model that mimicked the observed characteristics of these SSEs (Fig. 1a).The fault model represents the subducting plate interface as a single rectangular fault with a dip angle of 15 extending to 120 × 100 km in the strike and dip directions, respectively.The entire fault is divided into subfaults with sizes of 2 × 2 km; thus, the total number of subfaults is 3000.Each subfault was assumed to obey the following quasi-dynamic equation of motion (Rice 1993): where the slip velocity v i (t), shear stress τ i (t), and slip s i (t) (= dv i (t)/dt) are the temporally changing simulation variables at subfault i.The slip direction is fixed to the dip direction, which is roughly consistent with the relative plate motion of the Philippine Sea plate to the Amurian plate (Sella et al. 2002).The second term on the right side of Eq. ( 1) represents the radiation damping (Rice 1993) accounting for energy release owing to seismic wave radiation rather than the inertial term.The slip response function k ij represents a change in shear stress at subfault i due to a unit slip at subfault j and is calculated assuming a linear isotropic elastic homogeneous half-space (Okada 1992).The plate velocity v pl , shear modulus G, and shear wave velocity β are set as 6.5 cm/yr (Miyazaki and Heki 2001), 40 GPa, and 3.0 km/s, respectively.Note that this model implicitly assumes that the region outside the modeled rectangle fault moves at the plate velocity.
To express the slip behavior of SSEs, HN19 set a circular patch of SSE with a radius of R = 35 km, with its center corresponding to the center of the entire fault at a depth of 25 km.HN19 used the rate-and state-dependent friction law with the aging state evolution law described by Eqs. ( 2) and (3) (Dieterich 1979;Ruina 1983): where θ i (t) is a state variable representing the state of the sliding surface.τ 0 is the reference shear stress corresponding to the shear stress τ for steady sliding at a reference velocity v 0 ; however, these reference values were not explicitly set here because we only calculate the shear stress rate.The frictional parameters A i , B i , and (1) (2) L i that characterize the fault slip behavior at subfault i were assumed to be uniform inside and outside the circular patch.HN19 set the frictional parameters as (A-B, A, and L) = (− 50 kPa, 100 kPa, and 40 mm) and (70 kPa, 100 kPa, and 40 mm) for the subfaults inside and outside the SSE patch, respectively.These frictional parameters were determined based on the critical nucleation size R c : For an isolated patch, when R/R c < < 1, the calculated fault slip becomes stable, while an unstable slip occurs when R/R c > 1 for 0.5 < A/B < 1 (Chen and Lapsuta 2009;Hirahara and Nishikiori 2019).HN19 set the frictional parameters corresponding to R/R c ~ 0.47, such that the calculated slip behavior was similar to that of the longterm SSEs.The following calculations focus on only the frictional parameters inside the patch; the parameters outside the patch were set to the values used by HN19 throughout the study.
By combining Eqs. ( 1) to (3), the temporal evolution of the slip velocity v i (t) can be expressed as follows: By setting the initial conditions, we can numerically solve Eqs. ( 3) and ( 5) to obtain the temporal evolutions of two independent simulation variables, v i (t) and θ i (t), for each subfault using the fourth-order embedded Runge-Kutta method (Press et al. 1996).Consequently, HN19 successfully reproduced the observed characteristics, such as the recurrence intervals, durations, and maximum slip velocity, of the long-term SSE of the Bungo Channel.Hereafter, we focus on one cycle of the longterm SSE, and define the initial time as two years prior to SSE initiation, when dθ/dt was approximately zero.We fixed the initial values of the simulation variables at this initial time throughout the study.

GNSS observations
This study focuses on the Bungo Channel SSE that occurred from 2009 to 2011.Crustal deformations due to the SSE were monitored by the GNSS Earth Observation Network System (GEONET) operated by the Geospatial Information Authority of Japan.We used the daily GNSS time series from the 86 GEONET stations around the Bungo Channel that were used in HN19 after removing seven stations that were unavailable during the analysis period (Fig. 1a, Additional file 1: Table S1). (4) (5) The data period was from January 2006 to December 2011.The original time series were preprocessed using GipsyX-1.4software (Bertiger et al. 2020) under the precise point processing strategy with ambiguity resolution analysis.Station 0462 (Fukue) was used as a reference site (Fig. 1a).After preprocessing, we corrected offsets resulting from antenna replacement and large earthquakes by subtracting the differences between the coordinates, respectively, averaged for 10 days before and after the offset.We then rotated the two horizontal components to the trench-parallel (X-axis) and trench-perpendicular (Y-axis) components and subtracted the inter-SSE effect by fitting a linear function to the time series from 2006.5 to 2008.5.The standard deviation for each time series was calculated from the residuals of this linear trend fitting.Finally, we visually checked each time series to determine whether they included any ambiguous changes; if so, the station was not used in the subsequent analysis (see Additional file 1: Table S1).
Figure 1b shows examples of the GNSS time series, indicating that the transient signal was initiated in mid-2009 and continued for approximately 2 years.In the following analysis, we first assimilated the GNSS time series of the three components between 2008.5 and 2011.5 for fault slip monitoring.Subsequently, to test the predictability of the short-term evolution of fault slips, we changed the data period used for DA (hereafter referred to as the DA period) to 0.5, 1.0, 1.5, 2.0, and 2.5 years.In these trials, the initial time of the DA period was fixed to 2008.5, while the end was set to be 2009.0,2009.5, 2010.0,2010.5, or 2011.0 for 0.5, 1.0, 1.5, 2.0, and 2.5 years, respectively.

MCMC method
In this study, DA was conducted using MCMC, which is a technique used to obtain a realization or sample m from a target PDF p(m).In our problem, the target PDF is the posterior PDF of the model parameter p(m|d), where m consists of a set of three frictional parameters (A, B, and L) within the SSE patch, and d is the observation vector.Therefore, we aimed not only to optimize the frictional parameters, but also to evaluate their uncertainties.Using this information, the temporal evolution of fault slip can be predicted in a probabilistic manner.We adopted MCMC for the following reasons: MCMC is a simple and easy method to implement, as it solely requires forward models and does not employ any optimization technique such as the adjoint method (e.g., Kano et al. 2020) as long as computation time in a single forward calculation is short enough for realistic calculation; it is easy to modify the code when the forward model is updated; MCMC directly evaluates the posterior PDF of the model parameters of interest and provides probabilistic prediction, i.e., future evolution of slip velocity and the corresponding crustal deformation in our problem.However, this study mainly focuses on obtaining a single analysis of the model parameters that approximate the optimum posterior values rather than precisely evaluating the posterior PDF because the adopted simplified fault model provides a strong constraint on the evaluation of the posterior PDF, which we will discuss later.
Among the wide variety of MCMC methods, we used the Metropolis method, which is one of the wellknown versatile algorithms (Metropolis et al. 1953).The posterior PDF is calculated using Bayes' theorem p(m|d) = cp(m)p(d|m), where p(m) is a prior PDF of the model parameters, p(d|m) is a likelihood function, and c is a constant that is canceled out in the Metropolis method.Because we do not require constant c, we simply refer to the product of the likelihood function and prior PDF, p(m)p(d|m), as the posterior PDF and the realized value in a specific model parameter as the "posterior value" in the following results.It should be noted that in the fields of DA studies, we often define "cost functions" as metrics to be minimized.Since the minimization of the cost function is usually equivalent to the maximization of the posterior values, this study uses the posterior values of the probability as a metric from the view of the probabilistic approach.We assume two types of prior information: uniform distributions for all the frictional parameters, and those with constraints of R/R c < 0.78, that is, R c > 45 km, and 0.5 < A/B < 1.The latter prior PDF is introduced to avoid unstable slip and assumes that we know an SSE will occur in the target area.The likelihood function p(d|m) is defined as the product of the misfits between the calculated and observed crustal deformations, scaled by the observation errors with a normalization constant c': where d t includes the observed cumulative displacement on day t.s t contains the fault slip for all subfaults, and H is an observation matrix that converts simulation variables, or fault slip in this case, to observed quantities.Therefore, the product Hs t corresponds to the calculated cumulative displacement.In this study, the observation matrix H was calculated by assuming a linear isotropic elastic homogeneous half-space (Okada 1992).R t includes the observation errors described in the previous subsection, and N is the number of observation epochs.
The Metropolis method iteratively obtained the samples m from the posterior PDF p(m|d) using initial values of m 0 = (A, B, L) T = (100 kPa, 150 kPa, 40 mm) T , which are used in HN19.( 6) For k = 0, 1, …, K−1, the following steps were repeated: 1. Propose a candidate for the next sample m′ based on the following proposal distribution and the current sample m k : We assume Σ is a diagonal matrix.2. Determine whether the candidate m′ is accepted or not with an acceptance ratio of min( 1 To remove the effect of the initial values, the initial samples were not used.After this initial burn-in period, the resulting series of samples {m k } (k = 0, 1, …, K) emulates an objective sample set drawn from a target posterior PDF p(m|d).
For each iteration, we ran a forward simulation with the initial simulation variables and assigned frictional parameters to obtain the calculated crustal deformation, which was then detrended using the time series for the first two years.Subsequently, we compared the calculated displacement with the observations to obtain a sample.We first repeated these procedures for K = 10,000 iterations for the global parameter search with a standard deviation of 5.0 × 10 -1 kPa, 1.0 kPa, and 1.0 × 10 -1 mm for A-B, A, and L, respectively, in the diagonal components of Σ.Following this, we conducted an additional 10,000 iterations with values of 5.0 × 10 -3 kPa, 1.0 × 10 -2 kPa, and 1.0 × 10 -3 mm for A-B, A, and L, respectively, for the local search and obtained samples of the posterior PDF of the model parameters from the latter 10,000 samples.The number of iterations was manually determined based on the timing when the posterior values in the obtained samples show smaller change.

Fault sip monitoring results of the 2010 Bungo SSE
This section presents the results obtained by assimilating the GNSS data, including the entire SSE period from 2008.5 to 2011.5.Frictional parameters within the SSE patch were estimated as A-B = − 43.1 kPa, A = 77.5 kPa, and L = 45.9 mm (Table 1 and Fig. 2a-c).The uncertainty of each parameter was approximately six orders of magnitude smaller than the parameter value itself.This implies that a slight change in the parameters will greatly decrease the posterior PDF values, and the proposed model parameters will rarely be accepted in the MCMC (i.e., in step 2 in the previous subsection).The reason for the small uncertainty of the model parameters will be discussed in the next section.( 7) The calculated displacement time series using 100 randomly selected samples in Fig. 2a-c explained the observed GNSS data (Fig. 1b).Notably, all 100 calculated time series mostly overlap; thus, Fig. 1b indicates the small uncertainty of the calculated displacement, reflecting the small uncertainty of the parameters.Fault slips along the plate interface slowly accelerated in 2009 and lasted for ~ 2 years with a maximum slip rate of log(V/V pl ) ~ 0.55 (~ 11 cm/year) (Fig. 2d).These characteristics are roughly consistent with the kinematic inversion results of Yoshioka et al. (2015), demonstrating the reproducibility of the observed crustal deformation by DA.On the other hand, the calculated time series failed to capture the observed abrupt rise of the SSE in early 2010 (Fig. 1b and Additional file 1: Figure S1).It seems to be difficult to model this observed feature with high posterior values even if we change frictional parameters, suggesting a possible limitation of the assumed simple model, as we will show later in Fig. 4.

Short-term prediction of slip evolutions of ongoing SSE
We next investigated the short-term predictability of ongoing SSE by changing the DA period.This setting assumes that if we observe the currently ongoing transient signal, we can use DA to predict its future evolution over subsequent days to months based on the estimated frictional parameters.Fig. 3a, b and Table 1 summarize the DA and prediction results and the estimated frictional parameters.Fig. 3a, b, and Additional  − 116.7 (9.4 × 10 -4 ) 125.2 (8.9 × 10 -4 ) 39.3 (7.7 × 10 -5 ) 21.9 2.5 years a − 57.3 (3.9 × 10 -4 ) 67.4 (4.3 × 10 -4 ) 51.6 (1.8 × 10 -4 ) 61. file 1: Figure S1 indicate that all the results successfully explained the observed data during the DA period on average.However, the prediction results showed different behaviors depending on the DA period: when we assimilated data for a DA period of 2.5 years (orange line), the future evolution was predicted as a slow transient slip.However, when we assimilated data for DA periods shorter than 2.0 years (green, blue, and purple lines), the prediction resulted in a fast slip (Fig. 3b).This implies that without prior knowledge of observations showing the deceleration of a fault slip, or in other words, when we only have observations during slip acceleration, the results predict the occurrence of an earthquake.The estimated frictional parameters in the case of short DA periods correspond to unstable slip conditions of R/R c > 1 (Table 1).Following these results, we conducted similar DA trials with a prior constraint on the critical nucleation size R c to avoid unstable slip.The resulting time series calculated using the estimated frictional parameters (Table 1) predicted slow transient slip rather than unstable slip (green, blue, and purple lines in Fig. 3c).For cases with a DA period longer than 2.5 years, the proposed samples in the MCMC computation are not rejected by the prior constraints, and therefore, the results matched those without a prior constraint on R c .
Regardless of the prior constraint on R c , the uncertainties in the frictional parameters decreased with an increasing number of observation epochs (Table 1).This implies that the number of possible scenarios for future slip evolution decreased with data accumulation.To clarify this point, we analyzed possible scenarios for future slip evolution based on forward computations by varying the frictional parameters.We set a combination of frictional parameters A-B = − 210.0-30.0kPa, A = 10.0-150.0kPa, and L = 30.0-100.0 mm with increments of A-B = 10.0 kPa, A = 10.0 kPa, and L = 10.0 mm.Thus, we conducted 1862 forward computations, excluding the  and c (A-B, A). d The temporal evolution of the slip velocities along the cross-section is indicated by the dotted line in Fig. 1a sets of frictional parameters with which numerical integration could not be executed owing to strongly unstable behavior.In each computation, the posterior values were evaluated for all DA periods.Figure 4a summarizes all calculated time series or scenarios in each DA period in the Y-direction at station 1059.For trials with short DA periods, the time series with high posterior PDF predicted both fast and slow slips after the DA period.For those with long DA periods, only a few time series tended to predict fast slips (i.e., there is a reduction in the number of vertical red lines) and scenarios predicting SSEs were dominant.This characteristic can be clearly observed in Fig. 4b, which shows scaled posterior values as functions of the DA period on the horizontal axis, and in the final displacement in 2011.5 shown in Fig. 4a.The final displacement approximately corresponded to the mode of the calculated fault slips; the final displacement was large (> a few tens of cm) for fast slip, moderate (~ a few to 10 cm) for slow slip, and small (~ 0 cm) for stable slip.The scaled posterior values indicate that almost all the slip scenarios were possible for a DA period of 0.5 years.As the DA periods became longer, scenarios with fast slips were gradually rejected.Finally, those with high scaled posterior values favored a slow fault slip with final displacement of ~ 6 cm at station 1059, particularly in cases with DA periods longer than 2.5 years.
The uncertainty in the estimated frictional parameters was notably small relative to their values (Table 1), possibly due to the MCMC algorithm becoming trapped at a local maximum within the parameter space.It is generally difficult to judge whether the MCMC truly searches around the global maximum.Our result may search around the local maximum, although this local maximum is the best as long as within the searched range.Even if the global maximum better explains the data, our conclusion will be valid because our purpose is to find a model that yields to obtain predictions that better but not best match the observations.Additional file 1: Figure S2 shows the calculated time series using the sets of frictional parameters with the 30 largest posterior values among the 1,862 computations shown in Fig. 4a.This figure indicates that, in cases of DA period longer than 1.0 yr, the difference of posterior values log p (m|d) between the largest and the second largest is at least greater than 10, and thus the posterior values have one maximum at least within the searched range.We could conduct a number of forward simulations because the adopted forward model was simple and did not require a long computation time for a single calculation.However, this kind of discussion is generally unrealistic in the case of complex forward models, and therefore the implementation of parameter search algorithms reducing the possibility of being trapped into local optimum such as parallel tempering (Hukushima and Nemoto 1996), particle filter, and ensemble smoother with multiple data assimilation (Evensen et al. 2022), for example, would be a future update for a more global parameter search.Related to the local optimum problem, it is essential to confirm that the MCMC correctly samples the posterior PDF.Because we adopted the simplified fault model, this study mainly focuses on obtaining a single analysis of the model parameters that approximate the optimum posterior values rather than precisely evaluating the posterior PDF itself.For this reason, we simply determine the sample size based on changes of the posterior values.When we focus on a probabilistic prediction using a more realistic model in the future, metrics evaluating whether we obtain the sufficient samples to approximate the posterior PDF will be useful: the R-value (Gelman and Rubin 1992) is often introduced to quantitatively judge whether the samples have sufficiently converged; the effective sample size (e.g., Liu 2004) is a metric that objectively determines the sample size by considering the autocorrelations of the samples.Rather than the local maximum, the small uncertainties in the frictional parameters may be due to the strong constraint of the assumed simple numerical fault model, in which the location of the SSE patch is fixed, the frictional parameters are assumed to be spatially uniform within the patch, and the slip direction is fixed.These over-simplified assumptions inherently limit the representation ability in model outputs even if the model parameters are highly optimized through DA to partially mitigate these model errors.Future research should consider applying the method to a numerical model with high degrees-of-freedom, incorporating, for example, spatial heterogeneity in the frictional parameters, variable slip directions, and geometry of the subducting plate interface.In addition, we modeled SSEs using the rate-and state-dependent friction law with aging law as the evolution of state variables in Eqs. ( 2) and (3).It is pointed out that this mechanism seems inappropriate in explaining the observed features (e.g., Rubin 2008), and thus some other mechanisms such as dilatancy strengthening (Segall et al. 2010) and frictional behavior with cutoff velocity (Shibazaki and Iio 2003) were proposed to model SSEs.Integrating these mechanisms would be beneficial to examine the variability of prediction of spatio-temporal evolutions of fault slips based on the assumed mechanism in the future.It should be noted that MCMC is easy to implement with other models by substituting the forward equations.
This study assumed a uniform distribution of frictional parameters as a prior PDF.Bayesian estimation enabled the sequential update of the prior PDF, and the posterior PDF was eventually obtained after DA (Fig. 4b and Table 1).The posterior PDF could then be utilized as prior information when conducting DA for the next SSE.

Conclusion
We investigated the potential of DA for fault slip monitoring and the short-term prediction of SSE and applied the MCMC method to the 2010 long-term SSE in the Bungo Channel, Japan.The observed GNSS time series were quantitatively explained by the DA, independent of the DA period, even when using the simplified fault model proposed by HN19.For short-term prediction, information on the deceleration of slip velocities is necessary to appropriately constrain the evolution of the SSE; otherwise, the DA predicts fast fault slips as possible scenarios in terms of the posterior PDF.In this study, strong constraints were imposed on the fault model settings, assuming a planar fault geometry, the location of the SSE fault patch and spatial homogeneity of the frictional parameters.Even when using such a simplified model, the present results demonstrate the effectiveness of DA for monitoring and predicting the evolution of SSEs.The incorporation of a model with a high degreesof-freedom model could be promising for better fault slip monitoring and short-term prediction with high accuracy.

Fig. 1
Fig. 1 Target area and GNSS time series.a Top: the study area in southwest Japan is marked by the black rectangle in the lower image.The red square represents the Bungo Channel SSE fault region assumed by Hirahara and Nishikiori (2019).This study assumes uniform frictional parameters both inside and outside the SSE patch indicated in the red circle.The locations of GNSS stations are marked by green triangles, with the time series for stations A-D presented in b.The temporal slip evolution, denoted by the dotted line, is shown in Fig. 2d. a Bottom: the tectonic setting of the study area, in which the Philippine Sea Plate (PH) subducts beneath the Amurian Plate (AM) along the Nankai Trough.b Examples of GNSS time series at stations identified in a.The displacements in trench-parallel (X), trench-perpendicular (Y), and vertical (Z) components are represented by green, blue, and black dots, respectively.The red lines depict the calculated time series, derived from 100 frictional parameters sampled from the posterior PDF of the model parameters.Note that the red lines largely overlap, making it visually difficult to distinguish the 100 time series in each component.GNSS = Global Navigation Satellite System , p(m′|d)/p(m k |d)).If the candidate is accepted, we set m k+1 = m′; otherwise, m k+1 = m k .

Fig. 2
Fig. 2 Estimated frictional parameters and slip velocities.Frictional parameters were estimated by assimilating GNSS data between 2008.5 and 2011.5.Each panel shows 1000 samples obtained from the posterior PDF for a (A-B, L), b (A, L),and c (A-B, A). d The temporal evolution of the slip velocities along the cross-section is indicated by the dotted line in Fig.1a

Fig. 3
Fig. 3 Short-term prediction of crustal deformation.a Comparison of observed (circles) and calculated (colored lines) time series without prior constraints on the critical nucleation size R c in the trench-perpendicular component at station 1059.The calculated time series were computed using the frictional parameters with the highest posterior values by assimilating the observations indicated by the corresponding colored arrows.b An enlarged view of a in the vertical axis.c Same as a but with prior constraints on R c

Fig. 4
Fig. 4 Possible scenarios for future slip evolution.a Comparison of observed (black circles) and calculated (colored lines) time series, in the trench-perpendicular component at station 1059.These calculated time series are derived by varying a set of frictional parameters.The line colors correspond to posterior values calculated using the data in the DA period indicated by the black arrow in each subpanel.These posterior values are transformed to the form log10 (logp (m|d)) and then scaled to a range [0,1] within each subpanel.b Scaled posterior values as functions of the DA period in the horizontal axis and final displacement in 2011.5 shown in a in the vertical axis.The length of each horizontal black line corresponds to the scaled posterior values.Circles with warmer colors or longer horizontal lines correspond to scenarios with higher probabilities, while those with cooler colors or shorter lines indicate the scenarios with lower probabilities

Table 1
Estimated frictional parameters without and with prior constraints on the critical nucleation size Rc a Same values for both results (with and without prior constraints) Hirahara and Nishikiori (2019)ishikiori (2019)