Improving Detectability of Seafloor Deformation From Bottom Pressure Observations Using Numerical Ocean Models

We investigated ocean bottom pressure (OBP) observation data at six plate subduction zones around the Pacific Ocean. The six regions included the Hikurangi Trough, the Nankai Trough, the Japan Trench, the Aleutian Trench, the Cascadia Subduction Zone, and the Chile Trench. For the sake of improving the detectability of seafloor deformation using OBP observations, we used numerical ocean models to represent realistic oceanic variations, and subtracted them from the observed OBP data. The numerical ocean models included four ocean general circulation models (OGCMs) of HYCOM, GLORYS, ECCO2, and JCOPE2M, and a single-layer ocean model (SOM). The OGCMs are mainly driven by the wind forcing. The SOM is driven by wind and/or atmospheric pressure loading. The modeled OBP was subtracted from the observed OBP data, and root-mean-square (RMS) amplitudes of the residual OBP variations at a period of 3–90 days were evaluated by the respective regions and by the respective numerical ocean models. The OGCMs and SOM driven by wind alone (SOMw) contributed to 5–27% RMS reduction in the residual OBP. When SOM driven by atmospheric pressure alone (SOMp) was added to the modeled OBP, residual RMS amplitudes were additionally reduced by 2–15%. This indicates that the atmospheric pressure is necessary to explain substantial amounts of observed OBP variations at the period. The residual RMS amplitudes were 1.0–1.7 hPa when SOMp was added. The RMS reduction was relatively effective as 16–42% at the Hikurangi Trough, the Nankai Trough, and the Japan Trench. The residual RMS amplitudes were relatively small as 1.0–1.1 hPa at the Nankai Trough and the Chile Trench. These results were discussed with previous studies that had identified slow slips using OBP observations. We discussed on further accurate OBP modeling, and on improving detectability of seafloor deformation using OBP observation arrays.

(2) ΔP B (t) is pressure deviations around a reference pressure level according to the sea depth (P ref ), and is decomposed into geophysical and non-geophysical components. ΔP C (t) indicates the seafloor vertical crustal deformation. ΔP T (t) indicates the tides. ΔP O (t) indicates the oceanic variations. ΔP A (t) indicates the atmospheric pressure loading on sea surface. ΔP D (t) indicates the non-geophysical, instrumental drifts (Watts and Kontonyanis, 1990;Polster et al., 2009;Kajikawa and Kobata, 2019). ε(t) indicates residuals. The increase of 1 hPa in OBP is equivalent to the elevation of 1 cm in sea level or the subsidence of 1 cm in seafloor below OBP observation point.
When one investigates a particular phenomenon using OBP records, other components need to be estimated and removed. When dominant frequency of the target phenomenon is far from that of other noisy components, the target phenomenon is easily isolated using a Fourier analysis. However, time scales of slow seafloor deformation are comparable to those of oceanic variations. In particular, slow earthquakes with fault slips involve time scales of days to months and/or longer scales (Hirose and Obara, 2005). Vertical seafloor displacements associated with such slow slips are expected to be less than a few centimeters (Ohta et al., 2012) while oceanic variations are up to tens of centimeters in these time scales (Niiler et al., 1993;Donohue et al., 2010). Isolating slow seafloor deformation signals from OBP records is difficult.
Possible relationships between huge earthquakes and aseismic, slow earthquakes/slips at plate boundaries have been investigated and suggested (Ruiz et al., 2014;Tanaka et al., 2019;Baba et al., 2020). Seismological and geodetic efforts have been made to conduct campaign and continuous observations including OBP at the deep-sea bottom near the plate boundaries (Hirata et al., 2002;Hino et al., 2009;Kaneda et al., 2015;Wallace et al., 2016). It is necessary to develop suitable methods of separating slow seafloor deformation and oceanic variation in OBP observations. We may suppose that spatial scales of earthquake slips (<10 2 km) (Ohta et al., 2012;Sato et al., 2017) are mostly shorter than those of oceanic variations (>10 2 km) in time scales of days to months (Niiler et al., 1993;Donohue et al., 2010). When spatially dense OBP observations are carried out, based on this assumption, one can take approaches to cancel (unknown) oceanic variations by simply differentiating nearby OBP stations from a reference OBP station (Ito et al., 2013;Ariyoshi et al., 2014;Wallace et al., 2016;Sato et al., 2017), or by estimating spatially-common components using statistical and/or machine learning methods (Emery and Thomson, 2001;Hino et al., 2014;He et al., 2020).
On the other hand, removing oceanic and other components from single OBP record is a straightforward approach. Tidal signals have been precisely estimated by observed records or numerical predictions (Tamura et al., 1991;Egbert and Erofeeva, 2002). Lowerfrequency, non-tidal oceanic variations can be estimated by numerical ocean models, but their accuracy is not typically sufficient for detection of slow slips of less than a few centimeters in OBP observations (Fredrickson et al., 2019;Gomberg et al., 2019;Muramoto et al., 2019). Nevertheless, this approach of reducing ambient oceanic noises is fundamental, and is suitably applied for both multiple-station and single-station observations. During this decade, there have been extensive OBP observations installed for seismic/geodetic monitoring of shallow plate subduction zones around the Pacific Ocean ( Figure 1). Seafloor deformation due to slow slip has been found at a few specific regions and events (Ohta et al., 2012;Ito et al., 2013;Davis et al., 2015;Suzuki et al., 2016;Wallace et al., 2016). However, there was no report on general evaluations of detectability of slow seafloor deformation using OBP observations. Regional dependence of the detectability is not known. Accuracy of used numerical ocean models is not known in terms of OBP variations. In the present study, we investigate the current status of improving detectability of slow seafloor deformation using OBP observations. We utilize OBP data at several subduction zones around the Pacific Ocean, and available numerical ocean prediction models. The detectability is evaluated by root-mean-square (RMS) amplitude (i.e., standard deviation) of residual OBP time series at respective stations.

OCEAN BOTTOM PRESSURE OBSERVATIONS
We use OBP data listed in Table 1. Most of the OBP data are available at the Ocean Bottom Seismograph Instrument Pool (OBSIP) (http://www.obsip.org/). Data obtained at the Japan Trench were provided by Tohoku University (Hino et al., 2009;Hino et al., 2014). These were composed of campaign observations using pop-up-type autonomous gauges. DONET1/2 in Japan (Kaneda et al., 2015) and NEPTUNE-Canada (Barnes et al., 2015) are operated as seafloor cabled systems for long-term (>5 years) continuous observations. These observations are located at six subduction zones around the Pacific Ocean ( Figure 1; Table 1). The S-net cabled observatory was installed along the Japan Trench and the Kuril Trench, and the observed data have recently become available (Aoi et al., 2020). However, the S-net OBP data are currently undergoing quality control (Kubota et al., 2020), and are not used in the present study.
Far offshore OBP observations are more sensitive to offshore fault motions than onshore Global Navigation Satellite System observations. We investigate the OBP data obtained at sea depths greater than 500 m ( Figure 1). In addition, the utilized OBP data consist of approximately year-long continuous records, and we avoid using data with complicated drift or spike noise. In-situ OBP observations include various components as shown by Eq. 2. In the present study, oceanic variations at a period of 3-90 days  Bird (2003). Red and white circles denote used and unused stations, respectively. See Table 1 for detail. Isobaths of 500, 1,000, 2,000, and 4,000 m are shown.
are extracted using a band-pass filter to remove short-tide components and long-term instrumental drifts. The extracted oceanic variations are compared with numerical ocean models.

NUMERICAL OCEAN MODELS
Numerical ocean models with data assimilation have been developed to realistically predict and estimate oceanic states. Some of them have been available to end users. We use a single-layer ocean model and four multi-layer ocean models ( Table 2). The single-layer global ocean model (SOM) was developed by one of the authors Inazu and Saito, 2016). Multi-layer models are referred to as the ocean general circulation models (OGCMs). We use three global models (HYCOM, GLORYS, and ECCO2), and one regional model around Japan (JCOPE2M). Horizontal resolutions are 1/12°for SOM, HYCOM, GLORYS, and JCOPE2M, and ∼18 km for ECCO2. In terms of vertical resolution, these OGCMs contain 40-50 layers. Modeled oceans are driven by external forcing. The SOM is driven by wind and/or atmospheric pressure loading on the sea surface given by the JRA-55 reanalysis (Kobayashi et al., 2015). OGCMs are driven by wind forcing and heat/freshwater flux given by respective reanalysis datasets ( Table 2). Wind stress is the dominant driving force of ocean circulations with disturbances. The wind-driven ocean involves ∼10 2 cm in sea level variations (i.e., ∼10 2 hPa in OBP). Only the SOM takes into account the atmospheric pressure loading to drive the ocean.
The sea level mostly shows an isostatic response via gravity waves to atmospheric pressure loading (i.e., 1 cm elevation in sea level to 1 hPa depression in atmospheric pressure). This is known as the inverted barometer response (Wunsch and Stammer, 1997). When the inverted barometer response is established, the seafloor hardly feels the atmospheric pressure loading. The inverted barometer response is conventionally assumed especially in deep-sea regions. Heat and freshwater effects on the sea level changes are expected to be relatively small (<1 mm in sea level) in deep-sea regions (Ponte, 2006).

OBSERVED AND MODELED OCEAN BOTTOM PRESSURE
The OBP is shown as an integrated water pressure from sea bottom to sea surface, including atmospheric pressure loading above sea surface. The calculations of OBP are described below for SOM and OGCMs.
According to Inazu et al. (2012), SOM provides sea surface height (η) at each horizontal grid. OBP at a location (P SOM ) is simply given by the pressure due to sea surface height anomaly (η) plus atmospheric pressure loading on sea surface (P atm ): where ρ 0 , g, and H are constants of seawater density (1,035 kg/m 3 ), gravitational acceleration (9.81 m/s 2 ), and sea depth given at  each location, respectively. In the SOM, the wind and the air pressure at sea surface respectively drive the sea surface height and OBP. Figure 2 shows a scheme of sea surface height and OBP driven by wind and air pressure. The sea surface height driven by wind alone (η w ), which is directly projected to OBP (SOM w ρ 0 gη w ). The sea surface height (η p ) driven by air pressure (P atm ) shows mostly an inverted barometer response. However, there are certain amounts of departures from the inverted barometer in OBP (SOM p ρ 0 gη p + P atm ). Note that amplitude of SOM p is smaller than those of ρ 0 gη p and P atm , as indicated by Figure 2. The summation of OBP driven by respectively wind and air pressure (SOM w + SOM p ) is almost the same as the OBP simultaneously driven by wind and pressure (SOM w+p ), indicating that non-linear effects are negligible. Figure 3 shows spatial distributions of RMS amplitudes of OBP driven by wind alone (SOM w ), by air pressure alone (SOM p ), and by both simultaneously (SOM w+p ). SOM w+p is mostly comparable to SOM w . However, SOM p is also evident especially at the Southern Ocean and at marginal seas (Inazu et al., 2006;Ponte and Vinogradov, 2007).
The OGCMs provide sea surface height (η) and temperature/ salinity (T/S) at each vertical layer (depth: z), but currently do not include atmospheric pressure loading (i.e., P atm ). We calculate seawater density and pressure using these essential parameters. The EOS-80's equation of state for seawater (Fofonoff and Millard, 1983) is used to calculate the seawater density (ρ) from T, S, and hydrostatic pressure according to z. Vertical integration from sea bottom to sea surface is carried out to calculate OBP variations with a long-wave, hydrostatic approximation: Figure 4 shows an example of the calculation derived from HYCOM. It is evident that sea-surface height pressure (first term of Eq. 4) and sub-surface pressure (second term of Eq. 4) are inversely correlated. This compensation relationship between the sea-surface and the sub-surface pressure is typically predicted in a simplified two-layer ocean model (Gill, 1982;Unoki, 1993). The OBP is shown as an anomaly from the compensation.
As mentioned above, the calculated OBP from each ocean model is compared with observed OBP data at the period of 3-90 days. Figure 5 shows an example at the Hikurangi Trough. Good agreement is found between the observation and SOM driven by both wind and air pressure (SOM w+p ), as was reported by Muramoto et al. (2019). The agreement is considerably good at periods <10 1 days. A multi-layer model (GLORYS) shows a worse fit to the observation. GLORYS as well as other multi-layer models are driven mainly by the wind forcing, but not by air pressure. Here, we add an OBP component driven by the air pressure alone using SOM (SOM p ) to the GLORYS's prediction. The agreement between the observed and modeled OBP is notably improved especially at periods <∼10 days. This example indicates that the departure from an inverted barometer (SOM p ) is evident and should be incorporated to more accurately model the OBP at the period of 3-90 days even in the deep sea.

OBSERVED AND RESIDUAL OCEAN BOTTOM PRESSURE AT SIX SUBDUCTION ZONES
Comparisons between observed OBP and numerical ocean models are extended to the observation data around the Pacific Ocean ( Figure 1; Table 1). Results at representative stations at the six subduction zones are shown in Figures  6-11. Overall results in the respective regions are listed in Table 3. Detailed statistics are found in Supplementary Table S1.
The SOM and three multi-layer models (HYCOM, GLORYS, ECCO2) are evaluated for all the observations at the six subduction zones (Hikurangi Trough, Nankai Trough, Japan Trench, Aleutian Trench, Cascadia Subduction Zone, and Chile Trench). In addition to these four global models, JCOPE2M which covers the northwest Pacific Ocean is also tested for the observations around Japan (Nankai Trough and Japan Trench). The OBP data are first compared to that derived from the SOM driven by wind alone (SOM w ), and those derived from the multi-layer models (HYCOM, GLORYS, ECCO2, and JCOPE2M). Subsequently, OBP derived from SOM driven by air pressure alone (SOM p ) is added to the OBP derived from SOM w , HYCOM, GLORYS, ECCO2, and JCOPE2M, and those are subtracted from the observed OBP data. In this evaluation, we FIGURE 2 | Images of pressure variations at sea surface and sea bottom forced by wind and air pressure in SOM. ρ 0 gη w and ρ 0 gη p are sea level variations respectively forced by wind and atmospheric pressure (P atm ) in SOM. SOM w , SOM p , and SOM w+p are resultant ocean bottom pressure variations forced by wind, atmospheric pressure, and both simultaneously.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 598270 use the RMS reduction rate which is defined at the period of 3-90 days by: Figure 6 shows the result at the Hikurangi Trough in which the HOBITSS project  was carried out. Time series of eight stations from July 2014 to June 2015 were used in the present study. Mean RMS of these time series at the period of 3-90 days was 2.08 hPa. When SOM w , HYCOM, GLORYS, and ECCO2 were applied, the residual RMS amplitudes were 1.51, 1.66, 1.63, and 1.59 hPa, showing 20-27% RMS reduction rates by the respective models. When SOM p was added to these models (SOM w + SOM p , HYCOM + SOM p , GLORYS + SOM p , and ECCO2 + SOM p ), the residual RMS amplitudes were 1.21, 1.51, 1.30, and 1.46 hPa, indicating 27-42% reduction rates from 2.08 hPa. The contribution from SOM p was shown to be evident with 7-16% reduction rates. SOM w + SOM p provided the best RMS reduction of totally 42%, corresponding to a mean correlation coefficient of 0.82 between the observed and modeled OBP ( Table 3). Figure 7 shows the result at the Nankai Trough in which DONET1/2 (Kaneda et al., 2015) is operated. Time series of 45 stations from January 2018 to December 2018 were evaluated. Mean RMS of these time series was 1.19 hPa. When SOM w , HYCOM, GLORYS, ECCO2, and JCOPE2M were applied, the residual RMS amplitudes were 1.13, 1.63, 1.66, 1.29, and 1.42 hPa. This indicates that SOM w provided 5% RMS reduction, but other OGCMs hardly contributed to the RMS reduction. When SOM p was added to these models (SOM w + SOM p , HYCOM + SOM p , GLORYS + SOM p , ECCO2 + SOM p , and JCOPE2M + SOM p ), the   Only SOM w + SOM p provided the best RMS reduction of 17%, corresponding to a mean correlation coefficient of 0.58. SOM p contributed more to the RMS reduction than SOM w . Figure 8 shows the result at the Japan Trench in which Tohoku University carried out campaign OBP observations . Time series of 14 stations from May 2012 to November 2013 were evaluated. Mean RMS of these time series was 1.60 hPa. When SOM w , HYCOM, GLORYS, ECCO2, and JCOPE2M were applied, the residual RMS amplitudes were 1.49, 1.73, 1.57, 1.50, and 1.83 hPa. SOM w and ECCO2 provided 6-7% RMS reduction. Other OGCMs hardly contributed to the RMS reduction. When SOM p was added to these models (SOM w + SOM p , HYCOM + SOM p , GLORYS + SOM p , ECCO2 + SOM p , and JCOPE2M + SOM p ), the residual RMS were 1.35, 1.70, 1.43, 1.39, and 1.71 hPa. This shows that SOM w + SOM p and ECCO2 + SOM p provided RMS reduction of 13-16% (correlation coefficient of 0.54), indicating 7-9% contributions from SOM p . Matsumoto et al. (2006) reported that atmospheric pressure loading was useful to explain OBP variations at offshore northern Japan. A comparable result was confirmed using SOM with substantial number of data. Figure 9 shows the result at the Aleutian Trench in which the AASCE project (Barcheck et al., 2020) was carried out. Time series of 10 stations from June 2018 to March 2019 were used. Mean RMS of these time series was 2.02 hPa. When SOM w , HYCOM, GLORYS, and ECCO2 were applied, the residual RMS amplitudes were 1.84, 1.74, 1.81, and 1.78 hPa, indicating 9-14% RMS reduction by the respective models. When SOM p was added to these models (SOM w + SOM p , HYCOM + SOM p , GLORYS + SOM p , and ECCO2 + SOM p ), the residual RMS amplitudes were 1.78, 1.68, 1.76, and 1.75 hPa, indicating RMS reduction of 12-17% by the respective models. HYCOM + SOM p provided the best RMS reduction of 17% (correlation coefficient of 0.55). SOM p contributed to only 1-3% RMS reduction. Figure 10 shows the result at the Cascadia Subduction Zone in which the CI project (Toomey et al., 2014) was carried out, and NEPTUNE-Canada (Barnes et al., 2015) is operated. We used time series of 15 stations from October 2012 to June 2013 and from October 2014 to September 2015. Mean RMS of these time series was 1.93 hPa. When SOM w , HYCOM, GLORYS, and ECCO2, were applied, the residual RMS amplitudes were 1.80, 1.82, 1.72, and 1.77 hPa, indicating 6-11% RMS reduction by the respective models. When SOM p was added, the residual RMS amplitudes were 1.82, 1.86, 1.68, and 1.80 hPa, indicating that SOM p contributed to only less than 2% RMS reduction. GLORYS + SOM p provided the best RMS reduction of 13% (correlation coefficient of 0.51). Figure 11 shows the result at the Chile Trench in which the Chile PEPPER project (Tréhu et al., 2020) was carried out. Only four time series from May 2012 to March 2013 could be used. Mean RMS was 1.19 hPa. When SOM w , HYCOM, GLORYS, and ECCO2, were applied, the residual RMS amplitudes were 1.15, 1.35, 1.19, and 1.21 hPa. This indicates <6% RMS reduction by the respective models. When SOM p was added, the residual RMS amplitudes were 1.13, 1.38, 1.14, and 1.23 hPa. The contribution by SOM p to the RMS reduction was small (<3%). SOM w + SOM p and GLORYS + SOM p provided RMS reduction of 6-9% (correlation coefficient of 0.51).
The reduction of the residual RMS amplitudes provided by each ocean model with SOM p is summarized in Figure 12 and Table 3 (see also Supplementary Table S1). The RMS reduction rate of the time series is regarded as a proxy of model accuracy. As shown in Table 3, SOM w and/or OGCMs performed to represent observations at all the regions with RMS reduction of 5-27%. When SOM p was added to SOM w and OGCMs, the RMS reduction rates were up to 8-42% from the original oceanic variations at the period of 3-90 days. At the Hikurangi Trough, the Nankai Trough, and the Japan Trench, SOM w was relatively accurate (5-27% reduction) compared to other OGCMs, and SOM p was more effective (additional 9-15% reduction) than other regions. At the Aleutian Trench, HYCOM provided the best RMS reduction of 14%, and SOM p contributed to additional 3% reduction, with a total reduction of 17%. Other models provided 12-13% RMS reduction. At the Cascadia Subduction Zone, GLORYS provided the best RMS reduction of 11%, and SOM p contributed to additional 2% reduction, with a total reduction of 13%. At the Chile Trench, Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 598270 FIGURE 7 | Same as Figure 6, but for result at MRF22 at the Nankai Trough. Comparisons with JCOPE2M is added.
FIGURE 8 | Same as Figure 6, but for result at P03 at the Japan Trench. Comparisons with JCOPE2M is added.

Dobashi and Inazu
Bottom Pressure With Ocean Models FIGURE 9 | Same as Figure 6, but for result at LA34 at the Aleutian Trench.
FIGURE 10 | Same as Figure 6, but for result at M14 at the Cascadia Subduction Zone.
At the Hikurangi Trough, model representations were relatively good, and RMS reduction rates >20% were found by all the models (Figure 6). SOM was notably accurate as totally 42% RMS reduction including 15% contributions from SOM p . In the western Pacific regions (Hikurangi Trough, Nankai Trough, and Japan Trench), SOM performed better than other OGCMs to represent OBP. SOM employed the Japanese reanalysis data (JRA-55) as external forcing to drive the ocean, while other ocean models employed US or EU data ( Table 2). We speculate the accuracy of JRA-55 may be better than US/EU meteorological data in the western Pacific regions. In the eastern Pacific regions (Aleutian Trench, Cascadia Subduction Zone, and Chile Trench), HYCOM, GLORYS, and ECCO2 were slightly better than SOM or comparable to SOM. Differences of external forcing and/or incorporating vertical layers in OGCMs might contribute to these relative improvements.
As shown in Table 3; and Figure 7, OGCMs unexpectedly performed to increase the residual RMS amplitudes especially at the Nankai Trough. There is a strong western boundary current known as the Kuroshio current that frequently meanders in the vicinity of the Nankai Trough (Qiu and Miao, 2000;Ebuchi and Hanawa, 2003). Due to its meandering and associated mesoscale eddies, OBP also changes with several hPa changes in seasonal or longer time scales (Nagano et al., 2018). When ocean models wrongly represent Kuroshio variations, the residual RMS amplitudes possibly increase. In addition, if residual OBP time FIGURE 11 | Same as Figure 6, but for result at CP04 at the Chile Trench. series do not correlate well with SOM p , the subtraction of SOM p may also increase the residual RMS amplitude, which was found in some cases for the Nankai Trough (Supplementary Table S1). The SOM does not realistically capture the Kuroshio variations, but may weakly represent their features, potentially helping to reduce non-Kuroshio noises in OBP ( Figure 12).

DETECTABILITY OF SEAFLOOR DEFORMATION AND FURTHER ROOT-MEAN-SQUARE REDUCTION
The RMS amplitude of residual OBP time series is crucial in terms of detectability of the seafloor deformation. Based on the result of Figure 12, the evaluations are carried out for the residual RMS amplitude derived from numerical ocean models with and without SOM p . At the Nankai Trough and the Chile Trench, the residual RMS was smallest. When SOM p was added, the residual RMS amplitudes were 1.0-1.1 hPa. As shown in Figures 7 and 11, original RMS amplitudes were also small (1.1-1.2 hPa). The RMS reduction was 16% at the Nankai Trough, and 8% at the Chile Trench. The model accuracy was not so good at the Chile Trench, but the ambient oceanic noise was small there. At the Nankai Trough, possible slow slips were reported using the dense observatories of DONET (Suzuki et al., 2016). When one investigates slow seafloor deformation using OBP observations there, effects of the Kuroshio current should be carefully considered (Nagano et al., 2018). At the Chile Trench, there have only been a few relatively short observations (Tréhu et al., 2020). Future OBP observations at the Chile Trench may enable more detection of small fault slips compared to other regions.
At the Hikurangi Trough and the Japan Trench, the residual RMS amplitudes were both ∼1.5 hPa when SOM p was not considered. The residual RMS substantially decreased to 1.2-1.3 hPa when SOM p was added. Original RMS amplitudes were 1.6 hPa at the Japan Trench, and 2.1 hPa at the Hikurangi Trough. The model accuracy of SOM was best at the Hikurangi Trough, although the ambient oceanic components were larger than other regions. The expected detectability is perhaps equivalent between the Hikurangi Trough and the Japan Trench. Previous studies have used OBP observations to identify specific slow slips (Ohta et al., 2012;Ito et al., 2013;Wallace et al., 2016). Recent efforts to improve sensitivities to find transient/ramp steps in OBP records have been increasingly made using ocean models, statistical methods, and/or machine learning in these regions Gomberg et al., 2019;Muramoto et al., 2019;He et al., 2020).
At the Aleutian Trench and the Cascadia Subduction Zone, original RMS amplitudes were also large (1.9-2.0 hPa) compared to other regions. Contributions from SOM p to RMS reduction were small (<0.1 hPa). The residual RMS amplitudes in both the regions were ∼1.7 hPa, being still large when SOM p was added. Detection of small seafloor deformation may be difficult in these regions compared to other regions, due to larger ambient oceanic noises and lack of model accuracy. By using a local ocean state estimation based on HYCOM boundary condition, Fredrickson et al. (2019) carried out an effort to reduce RMS of OBP observations at the Cascadia Subduction Zone partially same as those used in our study, and they also had obtained results quantitatively similar to our study.
The inverted barometer response of sea level has often been assumed to be robust for periods greater than a few days, indicating SOM p has been often negligible. However, the results shown above suggested that SOM p was mostly non-negligible and provided 2-15% RMS reduction in OBP in time scales of days to months (3-90 days) at the six subduction zones. Note that mean air pressure at the sea level over the global ocean has seasonal amplitudes of less than a few hPa (Ponte, 1993;Wunsch and Stammer, 1997). Seasonal changes of the total freshwater volume in the ocean may contribute to <1 hPa changes in OBP over the global ocean . Sea level response to atmospheric Rossby-Haurwitz waves at a period of 5 days significantly deviates from the inverted barometer (Hirose et al., 2001;Mathers and Woodworth, 2004). Other mechanisms of departure from inverted barometer have been investigated as well (Stepanov and Hughes, 2006). The deviations from the inverted barometer (i.e., SOM p ) should be taken into account even in deep-sea regions when OBP observations are used to find small signals (<a few centimeters) of slow seafloor deformation.  Table 3. White numbers are correlation coefficients corresponding to the root-mean-square reduction rates. Gray and magenta indicate model's contribution without and with SOM p , respectively. At the Hikurangi Trough, the Nankai Trough, the Japan Trench, and the Chile Trench, SOM w + SOM p is shown as the representative. At the Aleutian Trench, and the Cascadia Subduction Zone, HYCOM + SOM p , and GLORYS + SOM p are respectively the representatives.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 598270 We utilized a single-layer ocean model of Inazu et al. (2012) to provide SOM p . Better representations of the oceanic variations induced by air pressure may be obtained by using OGCMs that represent realistic ocean interior. According to Ponte and Vinogradov (2007), air-pressure-induced OBP involves up to a 20% difference between stratified and non-stratified oceans. Comprehensive modeling including self-attraction and loading effects on seafloor (Stepanov and Hughes, 2004;Vinogradova et al., 2015) may help to improve representations of OBP variations. In addition to these reasonable ocean dynamics, boundary conditions of bathymetry and meteorological data should be also accurate to represent more realistic oceanic variations. Digital elevation models with global bathymetry are still going to be updated (Tozer et al., 2019). SOM employed the Japanese meteorological reanalysis data (JRA-55) which might be relatively accurate at Asia and western Pacific, compared to other regions. According to the result of Figure 12, substantial improvements at the Aleutian Trench, the Cascadia Subduction Zone, and the Chile Trench are expected if the accuracy of the meteorological data can be improved there. Gomberg et al. (2019) and Fredrikson et al. (2019) respectively utilized local OGCMs with HYCOM-based boundary conditions to reduce RMS amplitudes of OBP observations for slow slip detections. The RMS reduction might be more efficient if air pressure to drive the ocean is incorporated into their ocean models. Recently, Androsov et al. (2020) showed that atmospheric pressure loading has a role to improve correlation between observed and modeled OBP of both daily and monthly time scales at the Southern Ocean. When air-pressure induced components (eg, SOM p ) become widely available in addition to wind-driven OGCMs, such ocean models will become more useful for marine geophysics applications.
In the present study, we used a band-pass filter to compare extracted oceanic signals as well as to systematically remove longterm drift in OBP data at several subduction zones. However, band-pass filtering has potential to obscure a ramp change which may correspond to a slow seafloor deformation. When one tries to isolate slow seafloor deformation, band-pass filtering is not appropriate in most cases, and long-term drift should be carefully estimated and removed using other approaches. In this paper, we have focused on accurate estimation of oceanic variations. Note that if oceanic variations are better estimated especially at long periods, this will also help to better estimate long-term drift component. Drift estimation has been conventionally done by fitting a function (typically exponential + linear) to de-tided OBP time series (Eble et al., 1989;Polster et al., 2009). However, better drift estimation is expected by fitting the function to the OBP residual after removing both tides and longer-period non-tidal oceanic variations. Improved drift estimation will contribute to better detectability of seafloor deformation.

SUMMARY
OBP observation data at six subduction zones around the Pacific Ocean were investigated to improve detectability of slow seafloor deformation signals. Numerical ocean models were used to reduce RMS of oceanic components in OBP records at a period of 3-90 days. The ocean models included HYCOM, GLORYS, ECCO2, JCOPE2M, and SOM. The RMS reduction rate using Eq. 5 was used to measure accuracy of these models. The RMS reduction of 5-27% were provided using these models at the six regions. Departure from an inverted barometer response was calculated using SOM with atmospheric pressure loading, which was not considered in the currently available OGCMs. This component was effective with additional RMS reduction of 9-15% at the Hikurangi Trough, the Nankai Trough, and the Japan Trench. But the component was not so effective with only <3% RMS reduction in other regions (Table 3; Figure 12).
The RMS amplitudes of the residual OBP time series can be a proxy for detectability of slow seafloor deformation. The residual RMS amplitude depends on regions and accuracy of ocean models. As a result, the residual RMS amplitudes were 1.0-1.1 hPa at the Nankai Trough and the Chile Trench, 1.2-1.4 hPa at the Hikurangi Trough and the Japan Trench, and 1.7-1.8 hPa at the Aleutian Trench and the Cascadia Subduction Zone (Table 3; Figure 12). Detectability of slow seafloor deformation is expected to be better in regions with smaller residual RMS amplitudes.
These were the results for the data at specific observation periods. However, we expect that comparable features of the RMS reduction will be obtained for other data periods at the respective regions. Although there will be room to improve model accuracy, the present study evaluated the utility of currently available ocean models to reduce oceanic noise in OBP observations for detections of slow seafloor deformation in subduction zones. It is still not so easy to detect slow seafloor deformation of less than a few centimeters from single OBP station time series Muramoto et al., 2019). Further efforts will be suitably combined by differentiating time series and/or statistical method using a number of stations (Ito et al., 2013;Ariyoshi et al., 2014;Hino et al., 2014;Wallace et al., 2016).
Improving accuracy of numerical ocean models with data assimilation is essentially useful to accurately estimate and remove oceanic noises at each station. This study showed that OGCMs did not always perform better than SOM, and OBP changes induced by atmospheric pressure loading should be incorporated to enhance the representations of OBP variations at the period of 3-90 days. Improvements of accuracy of boundary conditions including external force (e.g, meteorological data) are also required as well as reasonable ocean dynamics for more precise estimation of OBP. Integrated approaches from oceanography and solid earth science are indispensable for improving future seafloor geodesy.

AUTHOR CONTRIBUTIONS
Both authors designed this research, carried out analyses, discussed the results, and wrote the paper.

FUNDING
This study was supported by JSPS KAKENHI Grant number 18K04654.

ACKNOWLEDGMENTS
We appreciate Prof. Ryota Hino, and Associate Prof. Yusaku Ohta of Tohoku University for kindly providing OBP data at the Japan Trench. Dr. Narumi Takahashi, Dr. Takeshi Nakamura, and Dr. Tatsuya Kubota of the National Research Institute for Earth Science and Disaster Resilience, Japan, kindly helped us to obtain OBP data of the DONET1/2 observatory. We thank Dr. Ruochao Zhang of the Japan Agency for Marine-Earth Science and Technology for providing comprehensive dataset of JCOPE2M. Discussion with Mr. Tomoya Muramoto of the National Institute of Advanced Industrial Science and Technology, Japan, was useful for analyzing OBP at the Hikurangi Trough. Comments from two reviewers and those from Dr. Laura Wallace, Guest Editor for the special issue, greatly improved the manuscript.