Strong‐Motion Broadband Displacements From Collocated Ocean‐Bottom Pressure Gauges and Seismometers

Dense and broad‐coverage ocean‐bottom observation networks enable us to obtain near‐fault displacement records associated with an offshore earthquake. However, simple integration of ocean‐bottom strong‐motion acceleration records leads to physically unrealistic displacement records. Here we propose a new method using a Kalman filter to estimate coseismic displacement waveforms using the collocated ocean‐bottom seismometers and pressure gauges. First, we evaluate our method using synthetic records and then apply it to an offshore Mw 6.0 event that generated a small tsunami. In both the synthetic and real cases, our method successfully estimates reasonable displacement waveforms. Additionally, we show that the computed waveforms improve the results of the finite fault modeling process. In other words, the proposed method will be useful for estimating the details of the rupture mechanism of offshore earthquakes as a complement to onshore observations.


Introduction
Recently, real-time ocean-bottom seismometer networks have been deployed at active tectonic margins for the dual purpose of both basic research and real-time monitoring (e.g., S-net and DONET in Japan, NEPTUNE in Canada, and OOI in the US; Aoi et al., 2020;Barnes & Team, 2007;Trowbridge et al., 2019).These networks enable us to obtain near-fault records associated with an offshore earthquake.When coseismic deformation is large enough, tsunamis occur due to such events and can damage coastal areas.Coseismic displacement records at near-fault stations are thus important because they have the potential to help in evaluating the earthquake source and its resulting tsunami quicklythis can improve the performance of early warning systems.In addition, they can be useful in estimating the details of the rupture mechanism of offshore earthquakes as a complement to onshore observations.The ocean-bottom record has complicated noise sources for seismic recordings of any kind (e.g., Hilmo & Wilcock, 2020;Webb, 1998), more so than in the onshore environment.
Ideally, we seek to obtain a displacement waveform through the double integrations of an acceleration or "strong motion" record.However, even in the simpler situation of onshore recordings, simple integration leads to unphysical results due to "baseline offsets" which are small shifts in the records (Iwan et al., 1985).Many schemes have been proposed to remove the drifts in onshore strong-motion records (e.g., Boore, 2001;Wu & Wu, 2007;Wang et al., 2011).
These can often succeed in baseline correction, but their convergence is not always guaranteed, and even more worrying, even if they do converge it is not always the case that they converge to the right final static offset (Boore et al., 2002).Additionally, baseline correction is difficult to apply in real-time settings.Alternatively, Bock et al. (2011) applied a Kalman filter approach to the records of the collocated high-rate Global Navigation Satellite System (HR-GNSS) and seismometer sites, and succeeded in estimating the displacement waveform.Lately, many papers expanded on the use of the Kalman filter approach with onshore records (e.g., Geng et al., 2013;Melgar et al., 2013;Niu & Xu, 2014;Tu et al., 2014;Zang et al., 2019).
In this study, we propose a new method to estimate the coseismic displacement waveforms from offshore strong-motion records.We base the approach on the general philosophy of Bock et al. (2011): using the Kalman filter and combining the strong-motion record with the collocated ocean-bottom pressure gauge (OBPG) records.However, we have added challenges in the offshore environment, since a coseismic OBPG record contains not only seafloor displacements but also seafloor accelerations, ocean acoustic waves, and tsunamis (Saito, 2013), we cannot use the same Kalman filter formulation employed for onshore records.
Section 2 explains the innovation of the Kalman filter approach proposed in this study.Then we apply it to synthetic data (Section 3) and real data for an M w 6.0 earthquake recorded offshore in the Nankai, Japan region (Section 4).Finally, in Section 5, the estimated displacement waveforms are evaluated through the finite fault modeling.
Note that our method focuses only on vertical displacement because OBPGs are only sensitive to physical phenomena that can be leveraged to calculate vertical displacement.Unless otherwise described all the variables in the later sections are the vertical component.

Kalman filter implementation
The Kalman filter is an optimal estimation method based on the state-space representation of a dynamic system (Kalman, 1960).It consists of the combination of two different models: (1) the dynamic model which captures the physics of the processes involved and describes the time evolution of the states, and (2) the observation model which establishes the relationship between measurements and the states.The models used in this study are expressed as: where  ,  ,  , and Ω are the displacement, velocity, tsunami height, and DC offset of acceleration.The variables of  are the properties of noise in estimation.The DC offset is defined by the difference in baselines between the observed acceleration   and the true acceleration   : Ω =   +   −   (Melgar et al., 2013).The vector [   Ω ]  is thus the state vector, the output or result of the Kalman filter estimator.In this study, we use the water-depth fluctuation ℎ =  −  and estimated tsunami height  ̃ as the observation (Eq.2).
Note that ℎ ̇ in Eq. 1 is the time derivative of ℎ; ℎ and  ̃ are independently estimated by any methods other than the Kalman filter as we'll discuss soon.The noise properties of each variable, , are considered to be the Gaussian with zero mean and a standard deviation, i.e.,  ~ (0, ), and independent of each other.
Because Eqs. 1 and 2 are for continuous time series, an extra step in the derivation is needed to discretize them (Lewis et al., 2008): where ∆ is the time interval of the dynamic model (0.01 sec in our case, which is the sample rate of the ocean bottom accelerometer);   and   are the Gaussian noise such as   ~ (0,   ) and   ~ (0,   ), respectively.The covariance matrices,   and   , are written as: Before applying the Kalman filter, we need to obtain ℎ and  ̃ independently of the ongoing estimation.In this study, they are estimated by the method proposed by Mizutani et al.
(2020), extracting the tsunami and displacement components from the coseismic OBPG records manuscript submitted to Geophysical Research Letters on a real-time basis, and the tsunami source model estimated by the tsunami waveform inversion using time-derivative waveform (Kubota et al., 2018), respectively (Fig. 1).When calculating  ̃, we assume that the tsunami occurs instantaneously, or, put another way that the rise time for the deformation of the seafloor and sea surface can be negligible.The tsunami inversion is  The covariance matrices,   and   , or the variances of each variable,   , are important tuning parameters of the Kalman filter.We propose them to be automatically calculated from the records by moving time windows.For the strong-motion records,    is the moving variance of  and   Ω is the absolute value of the moving average of ; the 5-sec window is used in both cases.
Since ℎ  is estimated as the difference between two filtered records with a moving taper window (Mizutani et al., 2020

Synthetic test
In this section, we confirm the validity of our method and assumptions with synthetic records.
First, the seafloor acceleration, velocity, displacement, and stress change were calculated by OpenSWPC, and then, GeoClaw calculated the tsunami caused by the seafloor movement under the non-linear long-wave approximation.
A two-layer velocity model was used for the synthetic test: a water layer with P wave speed of 1.5 km and density of 1.0 kg/m 3 ; a homogeneous sea-bottom layer with P and S wave speed of 7 and 4 km/s and density of 2.7 kg/m 3 .The thickness of each layer was 4 and 400 km, and the horizontal dimensions of the model domain were 200 × 200 km.The domain was divided into 0.5 and 1 km cells in the horizontal and vertical direction, respectively, and the time interval was 0.02 sec.The source was represented as a point source whose parameters were as follows: the moment magnitude was 7.0; the strike, dip, and rake were 0°, 45°, and 90° (a pure reverse fault); the depth was 10 km; and the rise time was 7.5 sec.We set this source at the center of the model domain (Fig. 2a).Fig. 2b shows the synthetic records directly above the source.The OBPG record had a large amplitude due to the ocean acoustic waves and sea-bottom acceleration.In acceleration records, we added the baseline offset based on the model of Iwan et al. (1985), i.e., strong ground motion causes an offset, and after that, a minor offset still remains.We defined the strong motion for the baseline shift as over 0.5 m/s 2 , following Iwan et al. (1985).The offsets during and after the strong motion were 0.07 and 0.0007 m/s 2 , which were visually adjusted.Although the baseline shifts slightly affected the acceleration record, they significantly altered the velocity and displacement records.Note that we applied the Kalman filter to the record at this station; others were used only for the tsunami inversion.
The estimated tsunami source model for  ̃ and the outcoming results of the Kalman filter are shown in Figs.2c and 2d, respectively.Although the input acceleration was contaminated by the noise, the displacement waveform estimated by our method agreed with the "true" one.
When calculating  ̃, our method assumes that a tsunami occurs instantaneously.To investigate the effect of the rise time, we conducted the same synthetic test with twice the rise time (15 sec).In this case, the threshold for the baseline shifts in the acceleration was defined as 0.1 m/s 2 .The resultant waveforms also agreed with the true one, particularly in the dynamic part of the displacement waveform (Fig. 2d).Since the seafloor and sea surface move simultaneously, it is difficult for the OBPG to observe any dynamic displacement components, that is, this manuscript submitted to Geophysical Research Letters agreement comes from the acceleration record.In other words, the method proposed in this study successfully combined the information of the collocated OBPG and strong-motion seismometer.

Application to real data
We used the records of Dense Oceanfloor Network system for Earthquakes and Tsunamis (DONET; Aoi et al., 2020), which is deployed off the coast of Kii Peninsula, Japan (Fig 3a) .
Each station of this network consists of an ocean-bottom seismometer and an OBPG.On 1 April 2016, an Mw 6.0 event occurred inside this network (e.g., Araki et al., 2017;Takemura et al., 2018).The OBPGs clearly observed the pressure change originated from the tsunami, and strong-motion accelerometer records observed significant shaking (the peak ground acceleration (PGA) was over 0.5 m/s 2 at near-fault stations) and contained clear baseline shifts (Kubota et al., 2018;Mizutani et al., 2020;Wallace et al., 2016).For the preprocessing of the acceleration data, we removed the pre-earthquake offset in the records by taking the mean of the record 10 sec before the earthquake.For the OBPG data, the ocean tide component and the offset were removed by the theoretical tide model (Matsumoto et al., 2000) and the mean of the 30 min of the pre-event record.When conducting the tsunami source inversion for  ̃, we calculated Green's functions with the GEBCO_2023 bathymetry (GEBCO, 2023) with grid intervals of 0.02° from unit sources which were set each 0.1°.Note that we excluded station KME18, the closest station to the source, from the tsunami source inversion as well as Kubota et al. (2018), although applied the Kalman filter to its record.The resultant model is shown in Fig. 3b.
Fig. 3c shows the displacement waveforms at stations KME18 and KME17, the second closest station to the source (the results of other stations are shown in Fig. S2).As with the synthetic test (Section 3), we succeeded in obtaining stable displacement waveforms.We, however, must pay attention to the residual displacement or DC component of the waveforms.
For example, the waveform at KME18 indicated a subsidence of about 10 cm.The same signal was observed in the OBPG record, which previous studies considered as a false signal (Kubota et al., 2018;Wallace et al., 2016).Since our method estimates the displacement via ℎ from OBPG records (Eq.4), the residual displacement was affected by such an unphysical offset in OBPG records.
To avoid this problem, we changed a method to estimate ℎ in Eq. 4. We now estimate ℎ from the same model for  ̃, and the variance   ℎ was also calculated in the same scheme as  ̃.
Note that this alternative method was applied only to the stations whose OBPG records might have been contaminated by unphysical offsets: KMA03, KME17, and KME18, as suggested by Kubota et al. (2018).The displacement waveforms obtained from this method were also stable (yellow lines in Fig 3c ) and the unphysical offsets could be removed.The residual displacements agreed with the one by the tsunami source inversion.

Finite fault estimation
To evaluate the utility of the displacement waveforms from Section 4, we next conducted a finite fault inversion and compared the resultant model with the one obtained from the inversion of tsunami records.We solved this linear inverse problem by the non-negative least square method (Lawson & Hanson, 1995) with spatial smoothing (Text S1 and Fig. S1).Green's functions were calculated by OpenSWPC and GeoClaw.In the seismic waveform calculation, we used the 3-D velocity structure model by Koketsu et al. (2012) with grid intervals of 0.25 km and a time step of 0.01 sec.In the tsunami case, we used the same model in Section 4, i.e., GEBCO_2023 bathymetry with 0.02º interval.
The number of subfaults used was 81, each with a subfault size of 5 × 5 km and a rise time of 3 sec.Here, we estimated only the slip amount at each subfault; the other fault parameters were fixed as the values of Wallace et al. ( 2016): the strike, dip, and rake were 215°, 5°, and 95°; the center of the fault model was at 33.385ºN and 136.434ºE, and at 9.8 km below the seafloor.We set the rupture speed to 2.1 km/s, 80 % of the S wave speed in this region (Kamei et al., 2012;Wallace et al., 2016).
Fig. 4a shows the result of the inversion using the displacement waveforms.To investigate estimation errors, we conducted a bootstrap method with 200 samples where we randomly selected stations at each iteration (Chernick, 2007).The M w calculated from this model was 5.9 and the VR was 68.7%.From the standard deviation of the bootstrap (right panel in Fig. 4a), a large slip patch beside the epicenter reflects the actual fault slip, while small one close to station KME20 is perhaps an estimation error.This result is consistent with the aftershock distribution detected by Japan Meteorological Agency (JMA) and the result from the tsunamionly inversion (the Mw and VR were 5.9 and 59%; Fig. 4b).We therefore conclude that the displacement waveforms estimated in Section 4 can be used reliably for studying offshore sources.The blue and orange lines represent the source time functions calculated by the models of (a) and (c).

Discussion
To combine displacement and tsunami data, we conducted an additional joint inversion (Text S1 and Fig. S1).The obtained fault model is shown in Fig. 4c.The VR were 57.1% and 44.1% for the displacement and tsunami; the estimated moment magnitude was 5.8.The model indicated that the slip propagated from the epicenter to the downdip direction.Compared to the model from the displacement waveforms only, the fault slip concentrated only around the epicenter.Its extent was also narrower than the model only using the tsunami.
Using the seismic records enables us to calculate the source time function (Fig. 4d).It indicated that a large rupture occurred at 3 sec after the origin time, and lasted 5 sec.The shorter duration compared to the one by the model using the displacement only reflects the slip distribution smaller than that.
Wallace et al. ( 2016) suggested that the aftershocks associated with this event occurred due to the afterslip immediately following the main shock because their fault model was separated from the aftershock cluster.Our model, on the other hand, agrees with the aftershock distribution and indicates that the aftershocks were caused directly by the main shock, that is, the afterslip may not be necessary for the aftershocks.
The aftershock distribution in Fig. 4 concentrates on the west of the slip.This is because JMA detected the earthquake location using only onshore stations.Araki et al. (2017) found slow-slip events after this earthquake located in the area between stations KME17 and KME18, different from the afterslip proposed by Wallace et al. ( 2016), which cover the north region of our fault model (Araki et al., 2017, Fig. 2b).In other words, our fault model can explain the aftershock distribution associated with this event sequence very well.

Conclusion
We proposed a new method to estimate the coseismic displacement waveform from collocated ocean-bottom strong-motion accelerometers and OBPGs.Through the synthetic test and the manuscript submitted to Geophysical Research Letters application to real data, we confirmed the displacement waveforms estimated by this method to be reliable.
On the other hand, at some stations close to the epicenter, the resultant waveform had a relatively large offset due to unphysical DC components in OBPG records.At present, we cannot remove such an offset automatically because it is difficult to model this offset, that is, cannot be included simply in the Kalman filter estimation.Although several studies investigated unphysical drifts in OBPG records, they focused on the static records (Chadwick et al., 2006;Hino et al., 2022).Clarifying the characteristics of coseismic OBPG records will improve the Kalman filter approach to ocean-bottom records.
The finite fault model that jointly inverted from both the displacement and tsunami waveforms showed improvements compared to the models estimated independently.In other words, the displacement waveform by our method can help us to reveal the details of the rupture process of offshore earthquakes.Figures S1 to S2

Introduction
Text S1 explains the details of the inversion methods in the main article.Fig. S1 shows the trade-off curves used to determine the weight in the inversion.Fig. S2 shows the displacement waveforms estimated by the proposed method at all the stations.
In the main text, we conducted three kinds of inversions: the tsunami source inversion (Section 2), the finite fault inversion (Section 5), and the joint inversion (Section 6).The first two solve the equation below: where , , , and  are the data vector, kernel matrix, spatial smoothing matrix, and model vector, respectively.The weight parameter  is determined based on the trade-off curve of the variance reduction (VR) and model variance.In this study, the variance reduction is defined as (Takemura et al., 2018): where  !"#$ () and  !$%& () are the observed and synthetic waveforms at station .
The trade-off curves of each result are shown in Fig. S1.
The tsunami source inversion is based on Kubota et al (2018).We first take the moving average with a time window of 60 sec and then apply a low-pass filter of 100 sec to the ocean-bottom pressure records.The time-derivative waveforms of them are used as the data and Green's functions.We set the record length to 25 min.The singular value decomposition is used to solve Eq.S1.
In the finite fault inversion, we solve Eq. 1 by the non-negative least squares method (Lawson & Hanson, 1995).For the tsunami data, we apply the same preprocessing as in the tsunami source inversion.For the displacement data, we apply a low-pass filter of 20 sec and use 30 sec from the origin time.
In the joint inversion, Eq. 1 is modified to:  . 2c, 3b, 4a, and 4b.The text at each point is  in Eq.S1, and the red circles represent the weight we used.Note that although the result of Fig. 4a comes from the bootstrap method, the trade-off curve is obtained by the single inversion.That is why the VR value is different from the main text.(b) Trade-off curves for the joint inversion (Fig. 4c).The left, center, and right panels are for   ,   , and  in Eq.S3, respectively.The text and color of each point indicates the weight and model variance.

Figure 1 .
Figure 1.The schematic flow of the proposed method.

Figure 2 .
Figure 2. (a) Sea-bottom residual displacement in the synthetic test.The Green star represents

Figure 3 .
Figure 3. (a) DONET stations used in this study.The blue, red, and green triangles are the

Figure 4 .
Figure 4. (a) Slips of the finite fault model (left) and their standard deviations (right) by the

Figure S1 .
Figure S1.(a)Trade-off curves used to determine the weight  in the inversions ofFigs.2c, 3b, 4a, and 4b.The text at each point is  in Eq.S1, and the red circles represent the weight we used.Note that although the result of Fig.4acomes from the bootstrap method, the trade-off curve is obtained by the single inversion.That is why the VR value is different from the main text.(b) Trade-off curves for the joint inversion (Fig.4c).The left, center, and right panels are for   ,   , and  in Eq.S3, respectively.The text and color of each point indicates the weight and model variance.

Figure S2 .
Figure S2.Same as Fig. 3c except that at all other stations.

Figure S2 .
Figure S2.(Continued) , Section 5.1),   ℎ is the sum of the variance of filtered records with the range of the taper's flat part (48 sec).
Kubota et al. (2018)rward are the kernel matrices used in the waveform inversion and forward calculation; [cov ] is the covariance matrix of the data used in the inversion, that is, the diagonal matrix whose elements are (   +  −1  )/∆, where    is the moving variance of the data with a time window of 60 sec.Note that we do not use    but (   +  −1  )/∆ because the waveform inversion ofKubota et al. (2018)uses time-derivative waveforms.