Extracting Near‐Field Seismograms From Ocean‐Bottom Pressure Gauge Inside the Focal Area: Application to the 2011 Mw 9.1 Tohoku‐Oki Earthquake

Recent studies have shown that ocean‐bottom pressure gauges (OBPs) can record seismic waves in addition to tsunamis and seafloor permanent displacements, even if they are installed inside the focal area where the signals are extremely large. We developed a method to extract dynamic ground motion waveforms from near‐field OBP data consisting of a complex mixture of various signals, based on an inversion analysis along with a theory of tsunami generation. We applied this method to the OBP data of the 2011 Tohoku‐Oki earthquake. We successfully extracted the low‐frequency vertical seismograms inside the focal area (f < ∼0.05 Hz), although those of the Mw ∼9 megathrust earthquake had never previously been reported. The seismograms suggested two dominant energy releases around the hypocenter. The seismic wave signals recorded by the near‐field OBP will be important not only to reveal earthquake ruptures and tsunami generation processes but also to conduct real‐time tsunami forecasts.

. However, near-field seismograms associated with the mainshock with a reasonable quality have not been reported. The high-sensitivity ocean-bottom seismometers (OBSs) installed off Miyagi  went off-scale and whole seismograms were not recorded. The strong motion accelerometers installed outside of the main rupture area (open triangles in Figure 1a) were dynamically rotated by the strong shaking (Nakamura & Hayashimoto, 2019). Although some near-source seismograms during past megathrust earthquakes have been recorded by onshore seismometers and GNSS, such as in the 2010 and 2014 Chile earthquakes (Madariaga et al., 2019;Vigny et al., 2011), the stations were located outside of the main rupture regions, where the permanent displacement was small.
During the mainshock, some ocean-bottom pressure gauges (OBPs) were installed around the main rupture area (inverted triangles in Figure 1). The deep-ocean OBPs often observe tsunamis, which have dominant frequencies lower than ∼0.01 Hz. Such tsunami data have been widely utilized, because tsunamis constrain the spatial extent of the seafloor vertical deformation (tsunami source) better than seismic waves (Kubota et al., 2018). This is attributed to tsunamis' much slower propagation velocity and there being a less significant tradeoff between the source dimension and rupture propagation velocity across the fault. Previous studies used the mainshock OBP data to investigate the mainshock tsunami generation process (e.g., Saito et al., 2011 [yellow contour lines in Figure 1a]; Baba et al., 2015;Dettmer et al., 2016;Gusman et al., 2012;Hossen et al., 2015;Maeda et al., 2011;Satake et al., 2013;Tsushima et al., 2011;Yamazaki et al., 2018). However, they did not utilize the OBPs installed inside the main tsunami source region where the seafloor uplift was extremely large (e.g., GJT3, Figure 1). This is mainly because there have been few near-field observation examples (e.g., Mikada et al., 2006) and the method to utilize the permanent deformation for tsunami modeling was not established. In this decade, the well-established method to utilize the permanent deformation for tsunami modeling was proposed  and many finite fault models using the OBPs inside the tsunami source have been obtained (e.g., Nemoto et al., 2019).
KUBOTA ET AL. Open triangles denote OBS stations by ERI. Black triangles are the F-net onshore seismometers. The white star is the mainshock epicenter  and the red CMT solution is taken from GCMT. Yellow contours denote the distributions of the initial tsunami height (Saito et al., 2011, 2 m interval). Gray dots and rectangular areas indicate the locations of the unit sources and the analytical area of the inversion analysis. The configuration of the unit sources in the space and time domains is schematically shown in the inset. (b) Pressure waveforms at GJT3. Gray, red, and blue traces are the de-tided, bandpass filtered (0.01-0.05 Hz), and lowpass filtered (0.01 Hz) waveforms, respectively. (c) Spectral amplitude before and after the mainshock (black and green, respectively), calculated based on Aki and Richards' (2002) definition. Passbands of the filters in Figure 1b  Our understanding of the ocean-bottom pressure change inside the focal area during tsunami generation has also progressed. The three-dimensional theory based on the tsunami generation mechanics showed the initial sea-surface distribution is not identical to the seafloor elevation (e.g., Kajiura, 1963;Saito & Furumura, 2009). Based on the three-dimensional theory, the pressure change inside the focal area includes a component related to vertically accelerating seafloor (Saito, 2013(Saito, , 2017(Saito, , 2019An et al., 2017;Filloux, 1982). This component was confirmed by the actual observations and numerical simulations (e.g., An et al., 2017;Filloux, 1982;Ito et al., 2020;Kubota et al., 2020;Kubota, Saito, et al., 2017;Matsumoto et al., 2012Matsumoto et al., , 2017Mizutani et al., 2020;Nosov & Kolesov, 2007;Saito, 2019;Saito & Tsushima, 2016;Webb, 1998). Additionally, the pressure change field includes the high-frequency ocean-acoustic waves related to the seawater elasticity (Lotto et al., 2019;Nosov & Kolesov, 2007;Saito & Tsushima, 2016). Recently the attempts to numerically simulate this complex pressure change field have started (e.g., Abrahams et al., 2020;Lotto & Dunham, 2015;Maeda et al., 2020;Saito et al., 2019). However, although the method to forwardly simulate more realistic OBP signals inside the focal has been established, it has also been reported that a simple bandpass filter cannot extract the seismic waves from the complex pressure change field (Saito & Tsushima, 2016). A method to appropriately decompose the OBP signal to the seismic and tsunami signals is not established yet.
The purpose of this study is to propose a method to appropriately extract the seafloor dynamic motion time series from the near-field OBP data inside the focal area. To achieve this, we attempt to decompose the OBP signals into seismic and tsunami wave signals based on a tsunami generation theory. Section 2 describes a theory of tsunami generation inside the focal area, the mainshock OBP data used in this study, and the procedure of our method. In Section 3, we show the results of the application of the method to the mainshock OBP data. Discussion and summary of this study are given in Sections 4 and 5, respectively.

Ocean-Bottom Pressure Inside the Focal Area
We represent the ocean-bottom pressure change inside the focal area as the sum of the contribution originating due to gravity (p gravity (t)) and that without gravity p non−gravity (t) (Saito, 2019): Supposing that the wave period is long, we may consider the seawater as an incompressible fluid. Also supposing that the sea-surface height change is small enough compared to the water depth and that the wavelength is much longer than the sea depth, p gravity (t) is approximately given by where ρ 0 = 1030 kg/m 3 and g 0 = 9.8 m/s 2 are the seawater density and gravity acceleration, and η(t) and u z (t) are the time series of the sea-surface height change (tsunami) and the seafloor vertically upward displacement, respectively. Hereinafter we refer to p hydrostatic (t) as the hydrostatic pressure change. The pressure change without gravity can be approximated as the dynamic pressure change, related to the action-reaction forces of the vertically accelerating seafloor, as (e.g., Saito, 2013;Saito & Kubota, 2019).
where h 0 is seawater depth. This relationship is basically valid at frequencies lower than the acoustic resonant frequency: f 0 = c 0 /4h 0 ∼ 0.05 Hz (c 0 : ocean-acoustic wave velocity). In this study, we attempt to extract the vertical acceleration d 2 u z /dt 2 from the pressure change p(t).

OBP Data
We use seven OBPs installed off Miyagi (green inverted triangles in Figures 1a and Hino et al., 2014) and two OBPs installed off Iwate (orange inverted triangles, Kanazawa & Hasegawa, 1997;Maeda et al., 2011). See Text S1 for the detail of these instruments. Station locations are listed in Table S1. We subtract the tidal components using the model of Matsumoto et al. (2000). We then apply a fourth-order Butterworth lowpass filter with a cutoff of 0.05 Hz in both forward and backward directions to reduce higher-frequency ocean-acoustic wave components. The cutoff of 0.05 Hz is determined considering the acoustic resonant frequency f 0 for the OBP at GJT3 (∼0.11 Hz). All records are resampled to 1 Hz after the filtering.
The de-tided waveform at GJT3 is shown in Figure 1b (gray trace). In Figure 1c, spectral amplitudes of the de-tided records before and after the mainshock are shown, which are calculated based on Aki and Richards' (2002) definition (time windows of 3276.8 s are used). High-frequency ocean-acoustic wave signals can be recognized even 1800 s after the origin time, and are dominant in frequencies higher than the acoustic resonant frequency f 0 ∼0.11 Hz. Dynamic pressure changes (Equation 3) are evident during the first few minutes, particularly for the frequency range 0.01-0.05 Hz (red traces in Figure 1b). Subsequently, low-frequency hydrostatic pressure changes (Equation 2) are also confirmed (<0.01 Hz, blue).

Extracting Ground Motions from OBP Data
This study attempts to extract the vertical acceleration d 2 u z /dt 2 in Equation 3 from the OBP data. In other words, our goal is to appropriately decompose the observed pressure change into its hydrostatic and dynamic components. To achieve this, we develop a method based on the inversion for the temporal evolution of the seafloor vertical deformation combined with the theory for ocean-bottom pressure inside the focal area described in Section 2.1. We represent the vertical displacement at the seafloor (u z (x,y,t)) by the superposition of functions U z,ij and  k, The basis function for the spatial distribution of the seafloor vertical displacement U z,ij (x,y) is given by which takes the maximum value at (x i ,y j ). The displacement time function    k t is given by where the function begins to increase at t = t k and reaches 1 after the duration T d . The coefficient m ijk in Equation 4 represents the displacement amplitude of the (i,j,k)-th function The hydrostatic pressure change at the n-th OBP located at (x n ,y n ) is given by n n n n z n n p x y t g x y t u x y t The first and second terms represent the pressure changes due to the tsunami and the vertical displacement at the seafloor at (x n ,y n ), respectively. The tsunami height η(x,y,t) is numerically calculated by solving the linear tsunami equation from the seafloor vertical displacement u z (x,y,t) (Equation 4). Since the p hydrostatic (x n ,y n ,t) is linear with respect to the seafloor displacement, we represent Equation 8 as the superposition using m ijk : We refer to   hydrostatic , , , ijk n G x y t as the hydrostatic pressure Green's function in this study, which is the hydrostatic pressure change at (x,y) excited by the unit vertical displacement of The dynamic pressure change at the n-th OBP located at (x n ,y n ) (Equation 3) is given by the displacement of Equation 4: By using Equations 8 and 9, we represent the pressure change at the n-th OBP excited by the vertical seafloor motions as where the Green's function G ijk (x,y,t) is given by which represents the pressure change at (x,y) excited by the unit vertical displacement of We estimate the displacement amplitude m ijk as model parameters in a linear inversion problem given by Equation 11, where the pressure change at the n-th OBP is used as the data. Using the estimated m ijk with Equations 8 and 9, the observed pressure change at the n-th OBP can be decomposed into the hydrostatic and dynamic components. The time history of the vertical acceleration can also be extracted using Equation 9, as KUBOTA ET AL.

10.1029/2020GL091664
In the same manner, vertical velocity and displacement can also be obtained. For example, vertical velocity is expressed as

Calculation of Green's Function and Inversion
We summarize the procedures of the Green's function calculation and the inversion analysis (see Text S2 for more detail). We suppose the x-and y-directions are along the trench-normal and trench-parallel directions, respectively. We distribute the spatial basis function U z,ij in an area of 220 × 270 km (gray dots in Figure 1a). We suppose the elliptical-shaped unit sources to be L x = 20 km and L y = 60 km and each of them overlaps with their adjacent ones (inset of Figure 1a). We also distribute the temporal basis functions  k with the temporal interval of  Δ 5 s t , during the 120 s from the origin time (inset of Figure 1a). Duration of  k is assumed as T d = 10 s. To calculate the hydrostatic Green's function, tsunami height is numerically simulated from the initial tsunami height distribution using the linear dispersive tsunami equation (e.g., Saito, 2019). We use the bathymetry data of GEBCO Bathymetric Compilation Group (2020), decimating to a spatial grid interval of 2 km. The input sea-surface height for the tsunami calculation is calculated from the unit seafloor displacement with the water wave theory assuming a constant depth of 6 km (Kajiura, 1963). The dynamic Green's functions are calculated, using the seawater depth h 0 for each station (Table S1). After the calculation of the Green's functions, the same filter as applied to the observation is also applied to the Green's functions.
In the inversion, we impose the constraints of the spatial smoothing (Baba et al., 2006) and spatial damping. The weights of each constraint are determined based on trial and error. The deformations are allowed to begin at t = 0 s. We use 3600-s time windows for the OBPs off Miyagi and 1800-s for the OBPs off Iwate.

Results
In Figure 2a, we compare the observed pressure changes at GJT3 with the synthesized ones (see Figure S1 for the other OBPs, Figures S2 and S3 for the spatial distribution of the total displacement and its temporal evolution). Figure 2b shows the extracted vertical accelerograms at the OBPs using the estimated model parameter m ijk in Equation 10. Compared to the observed pressure changes divided by ρ 0 h 0 (black traces), the extracted accelerograms (red traces) do not contain the low-frequency pressure signals due to the tsunami, which are evident after ∼120 s from the origin time. High-frequency pressure changes for the first 120 s are explained by the dynamic pressure components (green trace in Figure 2a) and the subsequent low-frequency pressure changes are modeled by the hydrostatic components (blue trace), and the overall pressure changes were explained very well by both pressure changes (red trace). From the amplitude spectra of the pressure change at GJT3 (Figure S4), we confirm that the calculated hydrostatic and dynamic pressure changes are dominant only in the low-and high-frequency ranges, respectively. In Figure 2b, we also plot the accelerograms of the onshore broadband strong-motion seismometer from the F-net (Okada et al., 2004, black triangles in Figure 1a) by gray traces. Although the arrivals of the main wave packet are delayed, the onshore seismograms are similar to the extracted ocean-bottom seismograms at the OBPs near each station (compare N.TYSF with TM1 and TM2,N.KSNF with P02 and P06,and N.KSKF with P03 and P07). We also show the vertical velocity and displacement waveforms in Figures 3a and 3b, respectively. The amounts of the calculated vertical displacements are surprisingly consistent with the observed pressure offset changes due to the permanent deformation (Figures 3c and S5 shows comparisons for longer time window). These comparisons indicate the validity of the extracted seismograms.
It is worth pointing out that the near-field seismograms inside the tsunami source where the vertical displacement was extremely large during the Tohoku-Oki earthquake had never been reported previously. In the accelerograms at the OBPs inside the main rupture area (GJT3, P08, and P09), two dominant positive pulses are confirmed (Figure 2b). The duration of the second pulse at GJT3 is relatively short compared to the first one, whereas the durations in both pulses at P08 and P09 are similar. From the velocity seismogram at GJT3, located ∼50 km landward from the trench axis, only one peak with a relatively long duration is confirmed (Figure 3a). On the other hand, at P08 and P09, located near the epicenter and ∼100 km from the trench axis, there are two velocity peaks at t ∼40 and ∼70 s (Figure 3a). These characteristics may reflect the rupture kinematics of the mainshock. One possible interpretation is that the rupture, or energy release, at the fault beneath P08 and P09, which are located near the epicenter, occurred twice. This feature is also suggested by the kinematic modeling of the mainshock from the regional or global seismograms (Lay, 2018).

Discussions
This study used a lowpass filter with a cutoff of 0.05 Hz to satisfy the condition that the seawater is considered as incompressible fluid. However, if the contribution of the seawater elasticity cannot be neglected in this frequency range, the extracted seafloor seismogram may be incorrect. To confirm validity of the extracted seismograms at f ≤ 0.05 Hz, we conduct a numerical simulation of the two-dimensional seismic wave propagation (Maeda et al., 2017, see Text S3 for the detailed settings). We assume the vertical cross-section passing through GJT3 along the trench-normal direction from the extended Japan Integrated Velocity Structure Model (Koketsu et al., 2012, top panel in Figure 4) and distribute point sources along the plate boundary. We apply lowpass filters with different cutoffs to compare the simulated pressure ( We adopted the inversion-based method to extract the ground motion signals from the OBP data. However, one might think that the bandpass filters are also capable of extracting the seismograms by removing the low-frequency hydrostatic components. In order to evaluate this, we investigate the accelerograms calculated based on a bandpass filter with passbands of 0.01-0.05 Hz, shown in Figure S6a (black dashed line traces). The bandpass filtered accelerograms seem to agree with those extracted by our approach. However, the waveforms do not agree at all when integrating to the displacement ( Figure S6b). This is because the permanent offsets are removed by the bandpass filter. Considering a slow rupture near a trench as in a tsunami earthquake (e.g., Lay et al., 2012), the megathrust earthquake rupture process possibly spans broadband frequency ranges. Because the spectral components of the mainshock ground motions possibly range into the low-frequency tsunami-dominant spectral bands, we must not use a highpass filter, which reduces the low-frequency components. It is essential to use a lowpass filter and to employ an inversion-based method with the tsunami generation theory to appropriately extract the broadband vertical ground motion including the low-frequency permanent offset component.
We could extract near-field seismograms from OBPs, which could never be achieved in the past when no OBP was installed inside the focal area and the tsunami generation theory was not established. By combining near-field OBP observation and the tsunami generation theory, it is expected that the parameters KUBOTA ET AL.   for the rupture kinematics and dynamics can be constrained more precisely, particularly for the subduction zone (e.g., Ide & Takeo, 1997;Kozdon & Dunham, 2014;Ma & Nie, 2019). In addition, developments in deep-ocean OBP observation enable us to capture the higher-frequency ocean-acoustic wave signals (Heidarzadeh & Gusman, 2018;Kubota et al., 2020;Webb & Nooner, 2016, Figure 1), which can be modeled by numerical simulation considering the seawater as the elastic body ( Figure 4, Maeda et al., 2017;Saito et al., 2019). Recent studies have started to simulate pressure changes inside the focal area based on the numerical simulations, accounting for both the hydrostatic and non-hydrostatic seismic components (e.g., Abrahams et al., 2020;Maeda et al., 2020;Saito et al., 2019). Besides, in deep-ocean measurements, it is still hard to control the installation environment and some studies have reported that the near-field OBS rotated due to strong shaking on the seafloor (Nakamura & Hayashimoto, 2019;Takagi et al., 2019). In such a situation, the near-field OBPs must produce powerful datasets to constrain the earthquake source information. Taking these facts into account, the high-frequency near-field OBP data should be more utilized to deepen our geophysical understanding of the subduction zone, as widely as the data from onshore and offshore seismic instruments.
Our approach utilizing dynamic pressure may also be applicable to practical real-time tsunami early warnings (e.g., Melger & Hayes, 2019; Tsushima et al., 2011;. Inside the focal area, the OBPs observe no hydrostatic pressure changes just after the origin time, because the sea-surface height change and seafloor vertical displacement are almost equivalent soon after the earthquake occurrence . If we utilize the dynamic pressure changes as vertical motion signals, which are dominant in the first few minutes, the accuracy of the tsunami forecast immediately after the earthquake rupture starts will be improved. In Figure S7, we apply the bandpass filter to the extracted seismograms (0.01-0.05 Hz, red traces) and compare to the seismograms expected from the observed pressure within that passband (black dashed KUBOTA ET AL.  traces). These two traces agree well with each other. This indicates the information of the band-limited velocity or displacement could be obtained from the bandpass-filtered pressure records, only in a few minutes from the focal time.

Conclusion
We developed a method to extract near-field seismograms from the OBP data inside the focal area. We applied the method to the near-field OBP data of the 2011 Tohoku-Oki earthquake to extract the ground motions inside the focal area, whereas the near-field seismograms during the Tohoku-Oki earthquake have never been reported yet. Our analysis successfully decomposed the OBP data into the dynamic pressure changes dominant in the first ∼120 s and the subsequent hydrostatic pressure changes due to tsunamis and permanent seafloor deformation. The extracted seismograms suggested that two dominant energy releases occurred beneath the OBPs near the epicenter. We confirmed the validity of the extracted seismograms based on the numerical seismic wave propagation simulation. Because the bandpass filter to reduce the low-frequency hydrostatic components also reduces the low-frequency ground motion components, our inversion-based method is essential to appropriately extract the ground motion waveform including the low-frequency permanent offset components. The high-frequency pressure change signals in the near-field OBP should be utilized more widely, for geophysical research as well as real-time tsunami forecasting.