Tsunami waveform inversion using Green’s functions with advection effects: application to the 2003 Tokachi–Oki earthquake

We explored nonlinear effects within the context of tsunami waveform inversion, wherein Green’s functions were linearly superimposed to estimate earthquake slips. We focused on these effects while developing a source model for the 2003 Tokachi–Oki earthquake off Hokkaido, Japan. A source model for this earthquake was developed based on linear tsunami waveform inversion using Green’s functions and tsunami waveforms observed at tide gauge stations. Subsequently, tsunami waveforms from the source were simulated at the stations using nonlinear long-wave theory and compared with those estimated by inversion. The comparisons demonstrated that the waveforms had a non-negligible discrepancy that was attributed to advection effects, even for the primary wave used in the inversion at the two stations. This result strongly suggests that advection effects should be considered in the source modeling of the 2003 earthquake based on tsunami waveforms observed by tide gauges. Based on these results, a new tsunami waveform inversion technique that incorporates linearly approximated advection effects and maintain the framework of linear tsunami waveform inversion using Green’s functions is proposed and applied. The proposed method successfully mimicked the advection effects during the 2003 tsunami, reproduced better tsunami waveforms, and developed a source model for the 2003 earthquake using these effects. The peak slip amount and seismic moment were greater in the source model with advection effects than those without the effects. This finding suggests that the values in the source models developed for other earthquake events without considering these effects may have been underestimated.


Introduction
Modeling and estimating the characteristics of past earthquakes are the most crucial areas of seismological studies and are equally important for estimating future earthquake characteristics.Inverse modeling, which utilizes observational data related to earthquakes, is a rational approach for estimating the characteristics of past earthquakes.Seismic waveforms, tsunami waveforms, and co-seismic crustal deformation data can serve as valuable inputs for inverse modeling (e.g., Robinson and Cheung 2010;Romano et al. 2010;Tanioka et al. 2004a).Tsunami waveform inversion methods have been applied to numerous tsunamigenic earthquake events worldwide, offering insights into their characteristics including location, dimensions, and rupture processes (e.g., Watada 2023;Tanioka and Yamanaka 2023).Analyzing tsunami waveforms, particularly those observed offshore, enables the estimation of the characteristics of the initial sea surface deformation (i.e., the initial tsunami) and the earthquake source strongly linked to these waveforms (e.g., Kubota et al. 2018).
The concept of linear tsunami waveform inversion (hereafter referred to as linear inversion) was introduced by Satake (1987).In this method, time-series data of sea surface changes (Green's functions) resulting from a unit slip in earthquake sub-faults were first simulated for stations, where a tsunami was observed as waveform data, based on linear theory.Subsequently, synthetic waveforms that best explained the tsunami waveforms observed at these stations were estimated using a linear combination of Green's functions.By synthesizing these waveforms, we can obtain the slip distribution of the earthquake as a solution to linear inversion, assuming a linear scaling relationship between the amount of slip and the magnitude of sea surface changes.In recent years, there has been a rapid increase in the number of tsunami observation systems in offshore areas, where nonlinear tsunami effects are negligible.Hence, utilizing offshore tsunami waveform data for linear inversion is a reasonable approach for estimating initial seasurface deformations or earthquake sources, because it guarantees the theoretical consistency of the solution obtained via linear combinations.
Tsunami waveforms in nearshore areas have been occasionally observed using tide gauges (tide gauge data) since the late 1800s (e.g., Kusumoto et al. 2020).Satake (1987) used tide gauge data rather than offshore observational data for linear inversion.When using tide gauge data for linear inversion, two factors require careful consideration: tide gauge response and nonlinear tsunami effects.Following Satake (1987), subsequent studies were conducted to investigate the effects of tide gauge responses and correct tsunami waveforms observed by tide gauges (e.g., Satake et al. 1988 and2010;Namegaya et al. 2009).Furthermore, tsunamis propagating over nearshore areas exhibit nonlinear effects unless they have a small amplitude.Consequently, the tide gauge data of a tsunami reflect both the tide gauge response and nonlinear effects.Although an approach that estimates a source considering the nonlinear effects in tsunami propagation was developed in a subsequent study (Pires and Miranda 2001), the linear inversion framework has been widely used for earthquake source estimation (e.g., Piatanesi and Lorito 2007;Satake et al. 2013;Gusman et al. 2015;Yamanaka and Tanioka 2021).In studies in which linear inversion was employed with tide gauge data, nonlinear effects were typically neglected or assumed to have an insignificant impact on the inversion results.The use of linear inversion with tide gauge data is contingent on the assumption that nonlinear effects can be disregarded.In this context, we raise two key questions regarding the utilization of linear inversion with tide gauge data: (1) When nonlinear effects are nonnegligible but not considered in linear inversion, what are the implications of this omission on the developed source model?(2) How can linear inversion be enhanced to accommodate nonlinear effects?
In light of the aforementioned challenges, we conducted a case study focusing on the 2003 Tokachi-Oki earthquake in Japan, which had a moment magnitude of 8.0 (Yamanaka and Kikuchi 2003;Tanioka et al. 2004a) and occurred along the Kuril Trench (Fig. 1).This earthquake caused a tsunami that affected the Hokkaido and Tohoku coasts.Tanioka et al. (2004b) promptly reported an overview of field survey results for the tsunami, and Tanioka et al. (2004c) published detailed survey results, including interviews, photos, and tide gauge data.Large earthquakes, such as those in 2003,1952, and the seventeenth century (e.g., Tanioka et al. 2004a;Hirata et al. 2003;Ioki and Tanioka 2016;Kobayashi et al. 2021), can be considered key elements for estimating the stress conditions within the southern part of the subduction zone along the Kuril Trench, and estimating their characteristics is essential for predicting future large earthquakes.Furthermore, large outer-rise earthquake events may occur after shallow megathrust earthquakes (e.g., Kodaira et al. 2021).As such, a careful investigation of the 2003 earthquake not only helps in understanding the rupture process but also in the prediction of potential large earthquakes in the region, including large outer-rise earthquake events.
Various data types, including global navigation satellite system (GNSS), strong-motion, tide gauge, and offshore pressure gauge data, were collected for the 2003 earthquake to estimate its characteristics, and several earthquake source models have been proposed in prior studies.For example, Yamanaka and Kikuchi (2003) employed teleseismic data, and Koketsu et al. (2003) used both GNSS and strong motion data in their inversion analyses to construct their source models.Tanioka et al. (2004a) conducted a linear inversion using tide gauge and offshore pressure data recorded during the tsunami.However, Tanioka et al. (2004a) assumed that nonlinear effects were negligible, even though the peak tsunami amplitude exceeded 2 m at the Tokachi tide gauge station.No studies have estimated the characteristics of the 2003 earthquake by accounting for nonlinear tsunami effects.To address the key questions presented earlier, we investigated the effects of the omission of nonlinear effects in linear inversion using tide gauge data on the proposed source model of the 2003 earthquake.In addition, our ultimate objective was to develop a tsunami waveform inversion method that incorporates nonlinear effects.

Linear inversion and nonlinear effects
Initially, a source model for the 2003 Tokachi-Oki earthquake was developed using linear inversion, following a method used by Tanioka et al. (2004a).The assumed source area in this study covered regions with significant slip, as estimated by Yamanaka and Kikuchi (2003) (Fig. 1).Tanioka et al. (2004a) estimated large slips within an assumed source area using sub-faults with lengths and widths of 40 km.In contrast, Tanioka et al.Yamanaka and Kikuchi (2003), and the area covered by the black rectangles outlines the assumed source area in this study (Tanioka et al. 2004a) (2004d) further subdivided the source area assumed by Tanioka et al. (2004a) using sub-faults with lengths and widths of 20 km for their inversion analysis.Kim et al. (2023) adopted a source setting similar to that of Tanioka et al. (2004d) and suggested that the 2003 earthquake might have been associated with larger slips than those estimated by Tanioka et al. (2004a and2004d).Tide gauge data, offshore pressure data, or both were included in these studies, along with their inversions.Based on these reviews, the tsunami data recorded during the 2003 event likely resolved the source characteristics with a resolution of up to 20 km.In this study, two scenarios were tested to assess the resolution capability of the inversion.We divided the assumed source area by eight sub-faults, each with a length and width of 40 km (Tanioka et al. 2004a) in the first scenario and 32 subfaults, each with a length and width of 20 km (Kim et al. 2023), in the second scenario.Each sub-fault within the assumed source area shared identical fault parameters of source location, dimension, and angles of strike, dip, and rake with those used by Tanioka et al. (2004a) and Kim et al. (2023).
We individually modeled the initial sea-surface deformations resulting from slips on the sub-faults as co-seismic vertical displacements using the method described by Okada (1992).The rupture velocity, duration, and rise time were not considered in any of the scenarios to estimate deformation (Kim et al. 2023).We combined three computational domains with spatial resolutions of 15, 5, and 5/3 arcsec within a nested grid system (Imamura et al. 2006).The domain with the highest resolution was constructed using data from the Geospatial Information Authority of Japan and the M7000 digital bathymetric chart provided by the Japan Hydrographic Association.The other domains were based on the General Bathymetric Chart of the Oceans (GEBCO), which originally had a resolution of 15 arcsec.We further modified the bathymetry in and around ports with tide gauge stations for better resolution, based on a marine chart published by the Japan Coast Guard.For the simulation, we used a linear long-wave (LLW) model based on Goto's (1991) model with linear and nondispersive assumptions: (1) where η denotes the change in water surface elevation from the mean sea level; M (= Uh) and N (= Vh) denote the fluxes in the latitudinal and longitudinal directions (λ and θ), respectively; U and V denote the velocities in the latitudinal and longitudinal directions, respectively; h denotes the still water depth; g denotes the acceleration due to gravity; R denotes the radius of the Earth; f denotes the Coriolis coefficient; t denotes time.The propagation of each initial sea surface deformation was then simulated using the LLW model and the resulting waveforms were stored as Green's functions (GFs).The simulation conditions are listed in Additional file 1: Table S1.This study used six tide gauge datasets from the Hanasaki, Akkeshi, Kushiro, Tokachi, Hachinohe, and Miyako tide gauge stations (Fig. 1).These data, collected at 1-min intervals, were identical to those used by Tanioka et al. (2004a).As part of a preliminary investigation, the GFs in the first and second scenarios were compared at the Tokachi station, which was closest to the significant slip area (Additional file 1: Fig. S1).Additional file 1: Fig. S1 shows that the variations in the GFs from various sub-faults were similar at the Tokachi station in the second scenario, and the GFs peaked at a similar time.In contrast, the GFs exhibited distinctive variations in the first scenario.These results indicate that the tide gauge data from the 2003 tsunami may not be capable of estimating the source characteristics with dimensions of 20 km.Consequently, we used GFs from subfaults with dimensions of 40 km.
As previously mentioned, tide gauge response effects may influence the observed tsunami waveforms.However, similar to Tanioka et al. (2004a) and (2004d), we did not account for these effects.Estimating the tide gauge response during the 2003 event was challenging, and we anticipated that the impacts were not significant at the stations based on the results of Satake et al. (1988).Using the observed waveforms and GFs, a linear inversion was performed based on the method of non-negative least squares (Lawson and Hanson 1995) to estimate the slip distribution on sub-faults with dimensions of 40 km (Fig. 2a).The time range used for the inversion was determined through trial-and-error (Fig. 2b). Figure 2a shows the estimated slip distribution with a peak slip of 5.2 m.Notably, Tanioka et al. (2004a), who conducted a similar inversion using tsunami data (tide gauge and offshore pressure data), estimated a peak slip of 4.3 m at the same sub-fault where we estimated the peak slip.This study and Tanioka et al. (2004a) applied different grid sizes to the finest computational domain.The resolutions used here and by Tanioka et al. (2004a) were 5/3 arcsec (approximately 50 m) and 4 arcsec (approximately 120 m), respectively.Thus, the difference in grid size contributed to the difference in peak slip values.Satake's (1995) finding that the use of small grids is effective in reproducing observed tsunami waveforms may also support this conclusion.
Subsequently, the propagation of a tsunami determined by the developed slip-distribution model was simulated based on a nonlinear long-wave (NLW) model.For the NLW simulation, we used the following governing equations, which were obtained by adding bottom friction terms to Goto's (1991) model: where P(= UD) and Q (= VD) denote the fluxes in the latitudinal and longitudinal directions, respectively; D (= η + h) denotes the total depth; n denotes Manning's roughness coefficient (0.030 and 0.025 m −1/3 s for land and water areas, respectively).A moving boundary condition was assigned at the boundary between sea and land to account for inundation.In the NLW simulation, the NLW model was applied only to the finest domain, whereas the LLW model was applied to the other domains (Additional file 1: Table S1).To investigate the nonlinear effects of several factors, we conducted three (4) simulations using a source model.Advection effects were taken into account in the three cases.However, in the first case, both bottom friction and inundation effects were neglected; in the second case, bottom friction effects were considered while inundation effects were neglected; and in the third case, both were considered.In cases where the inundation effects were neglected, a vertical wall condition was applied along the initial sea-land boundary (Yamanaka et al. 2023).
As shown in Fig. 2b and Additional file 1: Fig. S2, the tsunami waveforms generated by the LLW simulation (i.e., the synthesized waveforms) were closely aligned with the observed data during the primary waves.However, a moderate discrepancy was observed between the LLW and NLW simulation results for the primary waves at the Tokachi and Kushiro stations, whose data were employed in the inversion analysis.At these two stations, the initial peaks in the NLW simulations are smaller than that in the LLW simulation.These results indicate that for the 2003 event, nonlinear effects cannot be disregarded in the inversion analysis using tide gauge data.Neglecting these nonlinear effects may lead to underestimation of the slip amounts on the sub-faults.Notably, the waveforms obtained by the NLW simulation after the time used for inversion reasonably reproduced the observed waveforms better than those obtained by the LLW simulation (Fig. 2b and Additional file 1: Fig. S2).The inversion time range for the two stations covered the first peak time but did not reach the first zero-down crossing time.The discrepancy between the waveforms further increased with time after the end of the inversion.Therefore, these results suggest that the determined time range for inversion was reasonable for waveform reproduction using the LLW model.At the Hanasaki station, the peak tsunami amplitude during the primary wave was similar to that at the Kushiro station; however, the nonlinear effects were negligible.As such, the nonlinear effects were likely influenced by the characteristics of the tsunami and local topography.Additionally, bottom friction and inundation effects were small at all stations (Additional file 1: Fig. S2).At the Tokachi station, the NLW simulation without bottom friction and inundation produced the peak amplitude of 2.4 m during the primary wave, which was smaller than that of 2.9 m produced by the LLW simulation.The peak amplitude decreased from 2.4 m to 2.3 m owing to the bottom friction effects alone and further decreased to 2.1 m when both bottom friction and inundation effects were considered.Thus, we conclude that advection effects were most responsible for the discrepancy between the waveforms in the LLW and NLW simulations.

Modified linear long wave model
We aimed to incorporate the linear-approximated effects of advection terms into a linear inversion using the modified linear long-wave (mLLW) model developed by Yamanaka et al. (2023).The governing equations of the mLLW model in the spherical coordinate system can be introduced by adding new linear terms to the LLW model described as follows: where φ λ and φ θ are coefficients that vary in space.c denotes a calibration coefficient with a magnitude on the order of unity and is uniform in both space and time.
Referring to Yamanaka et al. (2023), we provide an overview of the characteristics of the mLLW model expressed in Eqs.(7-9).The second terms on the LHS of Eqs. ( 8) and ( 9) are linear terms that account for the (7 influence of the advection terms in the NLW simulation.
The values of advection terms during the primary wave are represented well by the parabolic curves of the flux values.By investigating the relationship between the values of fluxes and advection terms at each grid in the NLW simulation, the φ λ and φ θ values that reasonably approximate these parabolic curves are determined and multiplied with the fluxes.Furthermore, the φ λ and φ θ values are assumed to be consistently positive in attenuating fluxes in the mLLW simulation, leading to a decrease in the peak tsunami amplitude.However, in some regions, advection effects amplify the fluxes; in these regions, the values are set to zero.Avoiding negative φ λ or φ θ values enhance the stability of the mLLW simulation but may result in an underestimation of the amplification effects of the advection term in some grids.The influence of this underestimation of the amplification effects may also be accommodated by the calibration coefficient, c, which determines the relative magnitude of the dissipation effect of the advection term over the computational domain.The primary waves obtained in the NLW simulation are well reproduced by the mLLW simulation.However, the model remains linear.In contrast, in the NLW simulation, the peak time can be simulated earlier due to amplitude dispersion effects, which are nonlinear effects, compared to LLW and mLLW simulations.The waveforms estimated in the mLLW simulation are required to be time-shifted to match those estimated in the NLW simulation.
Using the NLW simulation results without bottom friction and inundation effects (presented in the previous section) with a sampling interval of 1 min and the same protocol as in Yamanaka et al. (2023), optimal values of φ λ and φ θ were determined in nearshore areas of the computational domains covering the Tokachi and Kushiro stations, where the nonlinear effects were not negligible (Fig. 3).Higher values of φ λ and φ θ indicate that the tsunami is significantly attenuated at the respective locations owing to the advection effects.Figure 3 shows similar results to those obtained by Yamanaka et al. (2023), wherein the magnitudes of φ λ and φ θ show significant horizontal variations across the areas and the highest values were observed at the mouth and edges of port breakwaters, where the magnitudes of the horizontal derivative of the flux appear to be large.These results affirm that in the NLW simulation, the 2003 tsunami was attenuated during its propagation across the areas, where the magnitudes of φ λ and φ θ were large owing to the advection effects.Additionally, it can be suggested that significant advection effects do not appear unless the computational topography effectively represented the port geometry.
Subsequently, the NLW model applied in the finest computational domain (Additional file 1: Table S1) was substituted with the mLLW model with optimal sets of φ λ and φ θ .The c value was set to 1.0 via trial-and-error.The GFs from the sub-fault sources were then estimated at the Tokachi and Kushiro stations based on these conditions and those estimated by the mLLW simulation were determined to be modified GFs (mGFs) (Additional file 1: Fig. S3).Moreover, the waveforms at the two stations were determined using the mGFs for the source model described in the previous section for comparison with those estimated using the NLW simulation (Additional file 1: Fig. S3).As shown in Additional file 1: Fig. S3 and Additional file 1: Table S2, the waveforms from the mLLW simulation were time-shifted by 49 s and 14 s for the Tokachi and Kushiro stations, respectively, to match those obtained via the NLW simulation within the gray time window area.The results shown in Additional file 1: Fig. S3 confirm that the waveforms estimated by the NLW and mLLW simulations match reasonably well.The waveforms estimated by the mLLW simulation, which considers advection effects, can be linearly superimposed, because the advection effects are approximated using linear terms (Yamanaka et al. 2023).Using this approach, a new source model of the 2003 earthquake reflecting advection effects is obtained by superimposing the mGFs at the two stations and the GFs at the other stations for each observed waveform.The linearity of mGFs enabled the use of the framework of linear inversion proposed by Satake (1987), even when considering advection effects, which was a unique characteristic of this approach.

Tsunami waveform inversion with advection effects
Slip distribution models with advection effects, denoted as S i = S i (E) (i = 0, 1, 2,…j), were sequentially developed in this study through an iterative process for slips, where S i (E) denoted the slip amount at a subfault (E).An outline of the inversion protocol is shown in Fig. 4, and the details are as follows.
The GFs(t, O) from unit sources was first estimated at each tide gauge station O = (λ obs , θ obs ) through LLW simulation.At i = 0, based on the linear inversion, a slip Optimal sets of φ λ (S i , L) and φ θ (S i, , L) were determined using the NLW simulation results, where L = (λ, θ) denotes location (Yamanaka et al. 2023).Assuming values for the calibration coefficient c, tsunami waveforms, η mLLW (t, O, S i ), are estimated through the mLLW simulation.In cases where the peak primary tsunami amplitude in η mLLW (t, O, S i ) successfully represents the amplitude of η NLW (t, O, S i ), the assumed value can be adopted as the value of c for the simulations in the iterative process.In this study, c was fixed at 1.0 in both computations for the Tokachi and Kushiro stations.The starting point of the iterative process was located here, and mGFs(t, O, φ λ , φ θ , c) from the sub-faults were then estimated by the mLLW simulation with the optimal sets.As outlined previously, the waveforms estimated in the mLLW simulation must be time-shifted to match those estimated in the NLW simulation, because the mLLW simulation does not account for amplitude dispersion effects.For correcting time delays in waveforms by the simulation, tsunami waveforms, η mLLW (t, O, S i ), were estimated by superimposing mGFs(t, O, φ λ , φ θ , c) for the source S i .After time-shifting η mLLW (t, O, S i ) to fit η NLW (t, O, S i ), the best time shift value, denoted as Δt (O), was determined at each station.In this study, the time-shift values had a resolution of 1 s (Additional file 1: Table S2).A new slip distribution model, S i+1 , can then be developed through superpositions of mGFs(t + Δt, O, φ λ , φ θ , c) at the Tokachi and Kushiro stations and GFs(t, O) at the other stations for η obs (t, O).Absolute values of differences in slips between the S i and S i+1 models, denoted as ΔS i , were calculated to determine if the estimated slips are convergent.We determined that the estimated slips were convergent if all values in ΔS i were smaller than 1% of the maximum slip in S i , and the S i+1 source model can be determined as the final source model with advection effects.Otherwise, we estimated the tsunami propagation from the source S i+1 by the NLW simulation, determined optimal sets of φ λ (S i+1 , L) and φ θ (S i+1, , L), and returned to the starting point of the iterative process with i counted up by one.
Figure 5a shows the developed source model with advection effects.In comparison between the source models with and without these effects, peak slips appeared on the same sub-fault, A2 (Figs. 5a and 2a).However, by considering the advection effects, the peak slip increased from 5.2 to 6.1 m (i.e., by approximately 17%).Furthermore, the slip on sub-fault A3 also increased from 2.5 to 3.5 m (i.e., by approximately 36%) (Table 1).In contrast, the slips at sub-faults B2 and B3 decreased slightly when accounting for advection effects.Consequently, assuming a rigidity of 40 or 65 GPa (Tanioka et al. 2004a;Kim et al. 2023), the seismic moment of the developed source model was 7.75 × 10 20 or 1.26 × 10 21 Nm (moment magnitude of 7.9 or 8.0), increasing from 6.67 × 10 20 or 1.08 × 10 21 Nm in the case without advection effects.Figure 5b compares the observed tsunami waveforms with those estimated by the NLW simulation using the developed source model with advection effects.The estimated waveforms reasonably match the observed waveforms.Notably, the underestimation of the amplitudes at the Tokachi and Kushiro stations by the NLW simulation, as shown in Fig. 2b, was improved (Fig. 5b).The observed peak amplitudes during the primary wave were 2.6 and 1.1 m at the Tokachi and Kushiro stations, respectively.Conversely, the peak amplitudes estimated from the source by the NLW simulation without considering the advection effects were 2.4 and 0.7 m, respectively, whereas those estimated with the effects reached Fig. 4 Outline of the inversion protocol 2.7 and 0.9 m, respectively (Figs. 2b and 5b).In addition, the slip distributions converged quickly in the iterative process.As shown in Table 1, three iterations were sufficient for the 2003 earthquake.These results demonstrate that our inversion method successfully developed a better source model by accounting for advection effects, indicating the efficiency of the proposed method.

Nearshore tsunami characteristics
Numerical simulations were conducted to validate the proposed source model.As presented in Sect."Linear inversion and nonlinear effects", the bottom friction and inundation effects were small for the source model developed using linear inversion (i.e., S 0 ).In this section, these effects are investigated again but for the final source model developed using mGFs (i.e., S 3 ).To this end, we conducted a simulation using the final source model for three cases.In the first case, both the bottom friction and inundation effects were neglected; in the second case, only the bottom friction effects were considered; in the third case, both of them were considered.
Figure 6 shows a comparison of the waveforms estimated by the NLW simulation for the three cases.As shown in the figure, the bottom friction and inundation had negligible effects on tsunami waveforms at most stations.However, at the Tokachi station, the peak amplitude during the primary wave decreased from 2.7 to 2.6 m owing to the bottom friction effects alone and further decreased to 2.4 m when both bottom friction and inundation effects were considered.The simulated peak amplitude of 2.4 m was smaller than the observed value (2.6 m) (Fig. 6) and comparable to that obtained by the NLW simulation (without bottom friction and inundation) using the source model without the advection effects (Fig. 2b).These results suggest that the slips  * The strike, dip, and rake angles were fixed at 230°, 20°, and 109°, respectively.The definition of each parameter was the same as that used by Tanioka et al. (2004a)

Source
Longitude (E) Latitude (N) Depth (km) Slip (m) in the developed source model were still underestimated, because the source estimation process did not account for inundation effects.Incorporating the inundation effects near the Tokachi station can be as important as considering the advection effects.
A plausible range of peak slips for the 2003 earthquake was investigated by comparing the simulated tsunami waveforms from the sources developed in this study with those of Kim et al. (2023).As previously mentioned, Kim et al. (2023) used a source geometry similar to that used by Tanioka et al. (2004d) and conducted a joint inversion using GNSS and offshore pressure data recorded during the 2003 event to estimate the slip distribution.The slip distribution model proposed by Kim et al. (2023) indicated a large peak slip of approximately 15 m in the vicinity of the Tokachi station.Figure 7a shows a comparison between the waveforms from the two sources based on the NLW simulation, accounting for bottom friction and inundation effects.At the Kushiro station, the tsunami generated by our source model was slightly underestimated compared to the observed waveform and smaller than that generated by the source model of Kim et al. (2023).Both models reproduced the observed waveforms at the station.However, the model developed by Kim et al. (2023) was likely better than our model in terms of reproducing the primary wave.In contrast, at the Tokachi station, the developed source model underestimated the observed tsunami, whereas Kim et al. (2023) significantly overestimated it.Figure 7b compares the inundation areas simulated using the two source models.Both models reasonably reproduced the inundation area inferred from the observations (Tanioka et al. 2004c); however, the inundation area produced by the source model of Kim et al. (2023) was greater than that of our source model.In the model of Kim et al. (2023), a significant slip area with values exceeding 10 m was concentrated near the Tokachi station.The slip was considered to have a substantial impact on the tsunami waveform at the station.Therefore, these results suggest that the large slip values reported by Kim et al. (2023), including the peak slip of 15 m, may have been overestimated.The large slip area of the Kim et al. (2023) model, with values exceeding 10 m, corresponds to subfault A2 in this study.To enable a comparison between the peak slips estimated in our study and those estimated by Kim et al. (2023), we calculated the average slip value of the Kim et al. 's (2023) model based on the area corresponding to subfault A2 in this study with a source dimension of 40 km.The calculated value was approximately 10 m, whereas the peak slip at sub-fault A2 in our study was approximately 6 m.Considering these results, along with the overestimation by Kim et al. (2023) and the underestimation in this study associated with inundation effects, the plausible range for the peak slip during the

Offshore tsunami characteristics
The ocean-bottom pressure gauges at KPG1, KPG2, TM1, and TM2 stations successfully recorded pressure changes associated with the 2003 earthquake and tsunami (Fig. 8).To investigate pressure changes due to the tsunami (hereafter referred to as the observed tsunami pressure), a low-pass filter with a cutoff period of 300 s was applied to the data after the astronomical tide components were removed using a tide prediction model (Matsumoto et al. 2000).In addition, a large-scale propagation simulation without a nested grid system was conducted based on the LLW model (Additional file 1: Table S1) to estimate the offshore tsunami characteristics for comparison with observational data.Figure 8 shows a comparison of the observed and estimated waveforms at the offshore stations.For comparison, the vertical axis in Fig. 8b was determined as the pressure head (H).With the hydrostatic pressure approximation, the pressure head can be calculated as H = (p − p 0 )/ρg = η − η 0 , where p denotes observed tsunami pressure, p 0 and η 0 denote observed tsunami pressure and change in water surface elevation at earthquake occurrence, respectively, and ρ denotes water density.As shown in the figure, the estimated waveforms were reasonably consistent with the observed waveforms.However, a moderate discrepancy was observed at KPG2.As mentioned previously, a source dimension of 40 km was used in this study, and it is possible that reasonable reproduction of the KPG2 waveform was challenging when this source dimension was used.
An important finding was that the waveforms from the source model without advection effects (i.e., S 0 ) were consistent with the observed waveforms (Fig. 8).This result indicates that the slip distribution model of S 0 closely resemble the actual slip distribution.The slip distribution patterns at S 0 and S 3 were similar (Table 1).These results suggest that the impact of nonlinear effects should have been considered for source modeling but was not large enough to drastically alter the slip distribution.In actual events, a tsunami with a large impact of nonlinear effects may occur (e.g., Yamanaka et al. 2023), and investigating whether the developed method is effective for such events will be future work.However, only moderate and small tsunamis can be fully observed using tide gauges, and nonlinear effects may not be substantial for such tsunamis.Therefore, the developed method is expected to be effective if it focuses on events in which a tsunami can be observed completely based on tide gauges.Note that it will be necessary to investigate whether waveforms determined by NLW simulation using a source model without considering advection effects reasonably reproduce waveforms observed by tide gauges to determine whether the source model is reasonable as an initial condition for the iterative process.

Conclusions
We conducted a case study of the 2003 Tokachi-Oki earthquake to investigate nonlinear tsunami effects and enhance the framework of linear tsunami waveform inversion originally proposed by Satake (1987) to account for nonlinear effects.The primary wave of the 2003 tsunami simulated using the NLW model was smaller than that simulated using the LLW model at the Tokachi and Kushiro tide gauge stations, and the results were primarily influenced by whether advection effects were considered.This result indicates that advection effects in the 2003 tsunami cannot be neglected and should be considered in source estimation.Considering these results, we proposed and applied a new tsunami waveform inversion method using modified Green's functions that incorporated the linearly approximated effects of tsunami advection as outlined by Yamanaka et al. (2023).The source was reconstructed using the proposed method and tide gauge data of the tsunami generated by the earthquake.In contrast to previous studies that used tide gauge data for tsunami waveform inversions without accounting for nonlinear effects (e.g., Tanioka et al. 2004a and2004d), we successfully incorporated advection effects into the source estimation for the 2003 earthquake.Consequently, the peak slip amount and magnitude of the developed source increased when the advection effects were considered.
As mentioned previously, many observation systems exist for tsunamis in offshore areas.However, the observational data from offshore tsunamis are limited to recent events.Although small tsunamis were occasionally observed by offshore pressure gauges before 2000 (e.g., Eguchi et al. 1998;Hino et al. 2001), the 2003 Tokachi-Oki event was the first major megathrust event for which near-field data from earthquake and tsunami source regions were recorded by seismometers and pressure gauges (Watanabe et al. 2004).Accordingly, many scholars have estimated the characteristics of large tsunamigenic earthquakes that occurred before 2000 based on linear tsunami waveform inversions using tide gauge data (e.g., Satake 1989;Hirata et al. 2003;Baba and Cummins 2005;Ioki and Tanioka 2017), with limited consideration of nonlinear effects.Considering our finding that linear inversion using tide gauge data may underestimate the peak slip and seismic moment, revisiting such a source may be necessary to improve the source model.
The inversion method proposed in this study can still be improved, because it does not consider other nonlinear effects, such as bottom friction and inundation.
In the case study of the 2003 event, these effects were not substantial.However, the inundation effects had a nonnegligible impact on the tsunami waveform simulated at the Tokachi tide gauge station.In future studies, the application of the proposed method to other events will be necessary to evaluate its performance, and the uncertainties of the source models developed in prior studies that did not consider nonlinear effects should be investigated.Furthermore, methods that consider the effects of bottom friction and inundation, in addition to advection effects, should be explored.

Fig. 1
Fig. 1 Source area of the 2003 Tokachi-Oki earthquake and locations of tide gauges (yellow triangles) used for tsunami waveform inversion.The contours represent a slip distribution model with an interval of 1 m estimated byYamanaka and Kikuchi (2003), and the area covered by the black rectangles outlines the assumed source area in this study(Tanioka et al. 2004a)

Fig. 2 a
Fig. 2 a Source model developed by the linear inversion; b comparisons between observed (black) and estimated waveforms using LLW (blue) and NLW (without both bottom friction and inundation effects) (red-dashed) simulations; the gray area indicates the time range used for the linear inversion

Fig. 3
Fig. 3 Spatial distributions of φ λ and φ θ in the nearshore areas of Tokachi and Kushiro.The yellow triangles indicate the tide gauge stations; the tick interval was 20 arc-s

Fig. 5 a
Fig. 5 a Source model developed with advection effects; b comparisons between the observed (black) and estimated waveforms determined via the NLW simulation (without both bottom friction and inundation effects) (red) and superposition of mGFs (blue-dashed); the gray area indicates the time range used for the linear inversion

Fig. 6
Fig. 6 Comparisons between observed (black) and estimated tsunami waveforms by the NLW simulation for three cases: without bottom friction and inundation effects (red), with only bottom friction effects (pink), and with both bottom friction and inundation effects (green).The gray area indicates the time range used for the linear inversion

Fig. 7
Fig. 7 Comparisons between tsunami characteristics estimated using the developed source model in this study and that in Kim et al. (2023) by the NLW simulation with bottom friction and inundation effects.a Tsunami waveforms at the tide gauge station; b inundation areas at the Tokachi area.The blue area and red point in panel (b) represent the simulated inundation area and tsunami observation point described by Tanioka et al. (2004c), respectively

Fig. 8
Fig. 8 Comparisons between estimated and observed waveforms at the offshore sites KPG1, KPG2, TM1, and TM2. a Locations of the tide gauges (yellow triangles) and offshore sites (red circles) and the source model with the advection effects; b waveforms at the offshore sites

Table 1
Fault parameters used and estimated in this study.