Precise Monitoring of Pore Pressure at Boreholes Around Nankai Trough Toward Early Detecting Crustal Deformation

In our recent study, we detected the pore pressure change due to the slow slip event (SSE) in March 2020 at the two borehole stations (C0002 and C0010), where the other borehole (C0006) close to the Nankai Trough seems not because of instrumental drift for the reference pressure on the seafloor to remove non-crustal deformation such as tidal and oceanic fluctuations. To overcome this problem, we use the seafloor pressure gauges of cabled network Dense Oceanfloor Network System for Earthquakes and Tsunamis (DONET) stations nearby boreholes instead of the reference by introducing time lag between them. We confirm that the time lag is explained from superposition of theoretical tide modes. By applying this method to the pore pressure during the SSE, we find pore pressure change at C0006 about 0.6 hPa. We also investigate the impact of seafloor pressure due to ocean fluctuation on the basis of ocean modeling, which suggests that the decrease of effective normal stress from the onset to the termination of the SSE is explained by Kuroshio meander and may promote updip slip migration, and that the increase of effective normal stress for the short-term ocean fluctuation may terminate the SSE as observed in the Hikurangi subduction zone.


INTRODUCTION
Inland dense networks of seismic and geodetic observations have revealed that slow earthquakes occur on the subduction plate boundaries in the shallower and deeper extensions of megathrust earthquake source regions worldwide (e.g., Schwartz and Rokosky, 2007;Obara and Kato, 2016). It is well known that slow earthquakes have longer duration times than regular earthquakes with the same moment release, and are classified into several types according to spatiotemporal scale (Ide et al., 2007). For instance, low-frequency tremor (LFT) and very low-frequency earthquake (VLFE) are detectable from seismometers because of their dominant frequencies of several hertz (Obara, 2002) and 10-100 s period (Ito et al., 2007), respectively. On the other hand, slow slip events (SSE) are thought to release aseismic slip for days to years (e.g., Ide et al., 2007), which would be detected from geodetic observations. Dense observation networks have revealed the relations between SSE, VLFE and LFT activity. Some of SSE accompany VLFE and LFT as "episodic tremor and slip" (ETS) (Dragert et al., 2001;Obara and Kato, 2016), where SSE fault zones overlap the source regions of VLFE and LFT.
For the shallower source region, slow earthquakes are close to the trench ( Figure 1). Recently, seismic observations on the seafloor such as at the DONET (Dense Oceanfloor Network system for Earthquakes and Tsunamis) (e.g., Kaneda et al., 2015;Kawaguchi et al., 2015) and S-net (Seafloor observation network for earthquakes and tsunamis along the Japan Trench) (Aoi et al., 2020) networks have been useful to monitor slow earthquakes in the shallower extension of subduction zones: e.g., LFT (Annoura et al., 2017;Tanaka et al., 2019) and VLFE (Nakano et al., 2018;Nishikawa et al., 2019). Even for the shallower extension far away from the inland observations, VLFEs are still detectable solely from inland observation networks (Takemura et al., 2019) because of their larger magnitudes than LFTs (e.g., Ide et al., 2007). On the other hand, SSEs are not detectable because they lack seismic signals and static crustal deformation decays with distance as r −3 , making their signals too small at inland GNSS networks far away from the trench.
Recently, the borehole observatories at C0002, C0010, and C0006 ( Figure 1) were successfully connected to DONET (e.g., Araki et al., 2017;Kinoshita et al., 2018). From March 2018, we can monitor crustal deformation at three boreholes along the Nankai Trough in real time. It has been well known that the pore pressure monitoring in boreholes is a useful tool to detect the crustal deformation driven by SSEs (Araki et al., 2017), where the inland observation networks and GNSS-A would not detect it because the fault slip is too small, estimated to be a few centimeters .  . Squares represent borehole stations. Labeled names (D-16, D-13, and C-11) are the nearby DONET stations for C0002 (blue), C0010 (red), and C0006 (yellow), respectively. Focal mechanism is obtained from Ariyoshi et al. (2021). (B) Seismic section showing the location of the boreholes C0002 (937-980 mbsf), C0010 (610 mbsf), and C0006 (453.5 mbsf), where mbsf meter below seafloor; VE vertical exaggeration. The average source location errors of sVLFEs are 5 and 2 km for horizontal and vertical component, respectively (Nakano et al., 2018). Spatial resolution of focal depth and fixed dip angle for the SSE is 1 km and 6°, respectively . (C) Borehole completion diagram at C0006. The solid magenta ellipses indicate the pressure gauges, where the top one (Pc0) is a hydrostatic reference pressure sensor used to remove oceanographic signals, BRT below rotary table, and MSL mean sea level. Figures (B,C) are modified after Figure 1 of Araki et al. (2017) and Figure 2 of Kinoshita et al. (2018), respectively.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 717696 Figure 2A shows the time series of pore pressure change at boreholes C0002, C0010, and C0006, shown in Figure 1, before and after the SSE in March 2020 (between the black dotted lines). In Figure 2A, the pore pressure change at C0002 is clearly seen as a ramp function. Noting that this characteristic is similar to that observed in the previous SSE in 2015 (Araki et al., 2017), we identified it as an SSE promptly after the occurrence of VLFEs. However, it took several months to estimate the fault model of the SSE .
This difficulty of deriving a fault model is due to the following circumstances. To avoid false recognition, an SSE is identified when significant signals appear at a minimum of two stations. For the SSE shown in Figure 2, we could not find significant signal at C0010 and C0006 due to low ratio of signal to noise in each. To overcome the difficulty, we extract the pore pressure change at C0010 due to crustal deformation by applying t-test between previous record without SSE and volumetric strain change driven by the SSE fault model .
Since SSE and VLFE activity is expected to become shorter recurrence interval (Matsuzawa et al., 2010) with greater moment release (Ariyoshi et al., 2014), it would be better to detect SSEs before the occurrence of VLFE toward early warning of megathrust earthquakes and tsunamis. From the analysis of previous SSEs observed from pore pressure in boreholes (Nakano et al., 2018;Ariyoshi et al., 2021), shallow VLFEs (sVLFEs) occur when the SSE migrates upward and its slip reaches the sVLFE source regions, and not so when the SSE region remains in the deeper part. However, we do not know possible factors of promoting slip migration upward. Since slow earthquake including SSE is thought to have stress drop much smaller than regular earthquake (e.g., Takagi et al., 2019), small perturbation such as oceanic oscillation may affect the generation process of SSE.
Recently, ocean modeling has been developed so as to conduct an ocean state nowcast/forecast system (JCOPE-T DA) that targets the coastal waters around Japan and assimilates daily remote sensing and in situ data (Miyazawa et al., 2021). JCOPE-T DA enables to extract the ocean fluctuation component driven by Kuroshio meander which impact on the seafloor pressure without tidal component under the baroclinic condition.
In this study, we develop the method of extracting crustal deformation from pore pressure in boreholes, and discuss the relationship between the occurrence of SSE and external perturbations such as atmospheric and oceanic pressure changes toward robust monitoring of plate coupling in subduction zones.

Removal of External Perturbation in Pore Pressure
There are three borehole observatories (C0002, C0010, and C0006) in the sedimentary wedge above the subduction plate boundary of the Nankai Trough ( Figure 1B). A schematic view of the station at C0006 is shown in Figure 1C. Pressure in the borehole (Pc1 to Pc3 in Figure 1C) is designed to be independent of bottom current by swellable packer. We have monitored the crustal deformation component of pore pressure changes by the sensors installed at hydraulically isolated depth, where fluid is confined by finely sediments with low permeability and water packer above the sensors. To monitor the crustal deformation Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 717696 close to the subduction plate boundary, we usually focus on the deepest pore pressure (Pc1) at each borehole. From previous studies, raw data of fluid pressure in boreholes is thought to also contain tidal and some external loadings other than crustal deformation (Davis et al., 2013) as shown in Figure 2B. To remove them from the pore-fluid records to obtain high precision pore pressure (Pp), we subtract pressure on the seafloor as reference (Pc0) with amplification correction factors, from pressure recorded beneath them in the borehole.
Focusing on Pc1 used in Figure 2, Pp at time t is described as where α is a damping correction factor. If there exits time lag (t 0 ) between Pc1 and Pc0, pore pressure at Pc1 is rewritten as where t 0 is expected as positive if the main loading is driven by oceanic fluctuations propagating from the seafloor downward. Previous studies adopt this strategy (e.g., Araki et al., 2017), but the detailed process of obtaining parameters of α and t 0 is not described. In the following sections, we investigate the spatiotemporal dependency of the parameters and their physical meaning.

Application of Nearby Dense Oceanfloor
Network System for Earthquakes and Tsunamis Sensors to the Seafloor Reference Figure 2 shows a drastic increase of the pore pressure at C0006 after May 2020 because of the drastic decrease of the pressure in Pc0, while there is no significant change at C0002 and C0010. During that period, VLFEs were not detected from broadband seismometers of DONET. If pore pressure is changed due to fault slip, temperature is expected to react due to fluid flux or frictional heat. However, seafloor temperature is quite stable (<0.02°C) at C0006 even during the increase of pore pressure after May 2020, where the laboratory experiment suggests that the temperature dependency factor on the pressure gauge is estimated as 3.1 hPa/°C (Machida et al., 2019). In addition, this increase has been going on now (over 1 year) at the same rate while temperature has been still stable. From these results, we conclude that the increase in pore pressure at C0006 is not due to crustal deformation.
If the abnormal increase of pore pressure at C0006 is due to the instrumental drift for the reference pressure gauge (Pc0) on the seafloor ( Figure 1C), it is reasonable to substitute nearby seafloor pressure gauges of DONET for the reference, where both pressure gauges of the boreholes and DONET use quartz crystal resonator sensors made by Paroscientific Co.. This means mechanical characteristics are expected to be similar to each other. In this study, we adopt seafloor pressure gauges at KMD16, KMD13, and KMC11 ( Figure 1A) for the reference of C0002, C0010, and C0006, respectively, toward real-time monitoring of crustal deformation even if the reference pressure gauges get out of order due to malfunction or maintenance operation. In this case, pore pressure is obtained as

Method of Extracting Crustal Deformation Component in Pore Pressure
On the basis of Removal of external perturbation in pore pressure and Application of nearby DONET sensors to the seafloor reference sections, we adopt the following process to estimate the values of α and t 0 for the pressure data of borehole and nearby DONET from the beginning of March to the end of June.
i) Synchronize the pressure data in the borehole and seafloor pressure of DONET by 1 s resampling. ii) Retrieve data segments of 708 h duration (the M2 tidal period). iii) Calculate the cross correlation function between Pc1_raw and the reference (Pc0_raw or P DONET _raw) on the seafloor. iv) Determine t 0 at the maximum value of the cross correlation coefficient. v) Calculate the value of α that minimizes the L1 or L2 norm of Pp between raw data of Pc1 and the reference (Pc0 or P DONET _raw) on the basis of Eq. 2 and Eq. 3. vi) Move the time window by one period of M2 tide (rough search) or one second (fine search). vii) Repeat the analysis starting at (ii) until the termination of the data is reached (end of June).
Since pressures are strongly affected by the tides, we adopt the time window for 57 intervals of the M2 tide, which is the almost 708 h and the least common multiple periods of M2 and S2 tides (59 intervals). If the time window is whole-number multiple of tidal period, we can avoid the dependency of time lag (t 0 ) on the length of time window for calculation of the cross correlation (Shirahata et al., 2019).

Ocean Model Product Used for Calculation of Atmosphere-Ocean-Induced Sea Floor Pressure Fluctuations
Practically, it is difficult to detide the observed seafloor pressure by using the theoretical tide within several hPa precision (e.g., . This is due to baroclinic condition and instrumental drift change excited by external factors such as seafloor current change (e.g., . To overcome this difficulty, we adopt the latest ocean modeling without tidal component. We calculate pressure fluctuations induced by atmospheric and oceanic variations at the seafloor using outputs provided by an ocean data assimilation system JCOPE-T DA (Miyazawa et al., 2021). JCOPE-T DA is derived not only from in-situ observations but also remote sensing data such as SST and SSH-A from satellite. Since the remote sensing covers wide range with high spatio-temporal resolution, the remote sensing data is powerful to enhance the reliability of JCOPE-T DA (Miyazawa et al., 2021). JCOPE-T DA is a data-assimilative ocean general circulation model developed from the Princeton Ocean Model with Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 717696 generalized sigma coordinate (Varlamov et al., 2015), covering 17°N-50°N and 117°E-150°E with horizontal 1/36 resolution and vertical 46 active levels. The model is driven by the atmospheric forcing such as sea level pressure, wind, atmospheric temperature and humidity provided from the National Center for Environmental Prediction Global Forecast System. The tidal forcing including the major eight constituents of K1, O1, P1, Q1, K2, M2, N2, and S2 is also applied on the model. The data assimilation is applied for correction of the model states based on various kinds of the observation data provided from the satellites, profiling floats, and ships. The hourly inputs of bottom pressure, temperature, salinity, and ocean current data are used for the analyses described in Oceanic fluctuation impact on the seafloor Section. The bottom pressure of JCOPE-T DA is calculated from the water density from sea surface to bottom and atmospheric pressure at sea surface. Temporal variation of bottom pressure variations calculated by JCOPE-T DA without external tidal forcing is represented as  Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 717696 where P a , ρ, η, η et , H denote atmospheric pressure, water density, sea surface height, tidal height obtained from the external tide model output (OTIS Oregon state university Tidal Inversion Software; Egbert and Erofeeva, 2002), bottom depth, respectively, and ρ 1 denotes the density of the first level grid. In the following sections, we investigate the time series of α and t 0 and apply the values that achieves the greatest noise reduction, and discuss the effect of external pressure changes by using JCOPE-T DA.

RESULTS
Temporal Change of α During the Slow Slip Event and Its Stability Before and After the Slow Slip Event Figure 3 shows the time series of factor (α) determined in each time window (of one period of M2 tide) at each station, for the borehole a) or DONET b) surface reference pressures. From Figure 3, the estimation of α at C0002 is generally stable, except for a sharp temporal increase during the SSE in both of the cases. This implies the pore pressure changed in the material surrounding the borehole during the SSE.
For C0010, the time function of α appears relatively stable before, during, and after the SSE. For C0006, where pore pressure change is expected to be negligible due to the SSE , it is quite stable in case of raw data of nearby DONET station and fluctuates in case of Pc0_raw for the reference seafloor pressure.
From these results, we conclude that the apparent increase of pore pressure at C0006 in Figure 2 is due to the instrumental drift at Pc0_raw, and that the instability as seen in Figure 2A at C0010 is not due to the malfunction of the reference pressure.
Relationship Between Temporal Change of t 0 and Tidal Perturbation Figure 4 shows the time series of time lag (t 0 ), where we apply the same operation to theoretical tide calculated by NAO.99Jb (Matsumoto et al., 2000) as shown in Figure 4B.  Figure 4B. As shown in Figure 4B and Figures 5G-I, time lag for the observation at nearby DONET station is greater than that for the theoretical one. This is because the pressure propagation of sea surface height change to the seafloor takes several seconds with decay due to baroclinic condition (e.g., Varlamov et al., 2015).
From Figure 4B and Figures 5D-F, time lag between Pc1_raw at C0006 and raw data of nearby DONET station KMC11 has a  Figure 5H) while the other is positive (Figures 5G,I). This is because KMD13 locates eastward from C0010 and is triggered tidal force earlier.
From Figure 4A and Figures 5A-C, the estimation of time lag between raw data of Pc1 and Pc0 is mostly zero for C0002 and C0010, whereas it fluctuates for C0006 because the operation of cross correlation does not work well due to the decrease drift for Pc0_raw ( Figure 2B). These results suggest that time lag between raw data of Pc1 and Pc0 is basically negligible for C0002 and C0010 while highly variable for C0006 because the drastic increase after May 2020 is larger than its tidal component.

Application of Fixed Parameters to the Entire Period of Data
In this section, we apply some patterns of α and t 0 fixed for the entire period of data. As described in Temporal change of α during the SSE and its stability before and after the SSE and Relationship between temporal change of t0 and tidal perturbation sections, the estimated values of α and t 0 at C0002 are affected by SSE and at all boreholes by the tidal component. From the time lag between raw data of Pc1 at C0006 and KMC11 in Figure 4B, there is time lag between observation (yellow) and theoretical tide (green) with periodic oscillations that are much longer than M2 tide component. This means that the assumption that tidal effects are negligible for the estimation of α and t 0 is incorrect and these parameters do not vary in accordance with the theoretical tidal calculation. Figure 6 shows time series of residual (L1 norm) and cross correlation coefficient with their scatter diagram. Figure 6C shows negatively correlation between residual and cross correlation coefficient. This means that the order of determination for α and t 0 does not affect the best solution. Figures 6A,B show residual and cross correlation became larger and smaller during the SSE. This implies that it is reasonable to find the fixed parameters so that time window does not include the duration time of SSE. From these results, we adopt the fixed parameters α and t 0 when the residual is the   Table 1. The others are the same as Figure 2.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 717696 lowest by the fine search for each reference setting. The obtained values are listed in Table 1. Here, we adopt L1 norm which can reduce the effect of spike change due to local earthquakes , but the results are the same as L2 norm except Pc1 for C0006. Figure 7 shows the time series of pore pressure with the fixed parameters in Table 1 in case that the reference pressure on the seafloor is assumed as Pc0 and P DONET _raw. From Figure 7B, we find that the apparent increase of pore pressure at C0006 is removed by substituting Pc0_raw to KMC11 as the reference, while time lag between raw data of Pc1 and Pc0 at C0006 cannot cancel out the tidal component in Figure 7A. For the overall of C0002 and C0010, the pore pressure aligned by Pc0_raw is similar to that aligned by raw data of nearby DONET stations. Figure 8 shows the amplitude spectrum of the pore pressure at Pc1 for each borehole station, focusing on M2 (12.42 h) and S2 (12 h) tides. For C0002, tidal component of M2 is removed by 30-40% for both Pc0 and KMD16 while S2 is increased by 20% for Pc0. These results suggest that the combination of M2 and S2 are totally decreased about 10-20%. For C0010, tidal component is removed by 10-30% for both of Pc0 and KMD13. Considering that the amplitude around M2 and S2 tides (0.52 and 0.5 days) is a significantly large peak, we think that the fluctuation at C0010 in Figure 2 and Figure 7 is excited by tidal components with nonlinear time lag possibly due to failure of packing water in Figure 1C. For C0006, tidal component is removed by 60-90% for KMC11.
Toward monitoring of pore pressure in real time, we determine the fixed value of α and t 0 under the present condition. From the values of minimum residual and maximum cross correlation coefficient in Table 1, the reference pressure at C0002 and C0006 should be Pc0 and KMC11, respectively. For C0010, the values of the residual and cross correlation coefficient seem almost the same between Pc0 and KMD13. In addition, both of time lag (t 0 −15 and −60 s) is out of mode, average and median values as shown in Figure 5E. This indicates that the parameter values at C0010 in Table 1 might be determined by chance. To investigate the robustness, we also plot the amplitude spectrum in case of Pc0 without time lag (t 0 0) in Figures 8cd, which shows that there seems no significant difference from the case with time lag. Therefore, it might be practical not to apply time lag in case of Pc0 at C0010, which is the same as the best condition for C0002.

Crustal Deformation Process Close to the Nankai Trough
Applying KMC11 to the reference seafloor pressure, we get the pore pressure at C0006 with high quality as demonstrated by Figure 7B and Table 1. Using this result, we discuss the crustal deformation there more detail. Figure 9 shows the close up of pore pressure history at C0006 before and after the SSE occurred in March 2020, which shows there exists compression about 0.6 hPa during the SSE. Considering that the onset and termination of the pore pressure change at C0006 is consistent with the duration time of the SSE determined from the pore pressure change at C0002 , we assume that this increase is due to crustal deformation driven by the SSE. Giving the conversion factor of C0006 as 5.0-6.0 kPa/μstrain (Davis et al., 2013;Araki et al., 2017;Ariyoshi et al., 2021), we roughly estimate the volumetric strain as 0.01 μstrain. This amount is about 10 times greater than expectation (0.0012 μstrain) from the fault model derived by Ariyoshi et al. (2021).
Considering the plate bending near the Nankai Trough as shown in Figure 1B and obtained dip angles of focal mechanism for the sVLFE in Figure 1A are in the range of 0-6°, we think that this discrepancy of the volumetric strain amount is mainly due to lower dip angle than the fault FIGURE 9 | Close up of pore pressure history at C0006 before and after the SSE in March 2020 by applying KMC11 to the reference seafloor pressure. Vertical double arrows represent change amount during the SSE.  Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 717696 models (A-1 and A-2+ in Figure 1B) estimated by Ariyoshi et al. (2021), which assumes a planar fault dipping constantly at 6°. To roughly estimate the valid dip angles for A-1 and A-2+ models, we calculate the volumetric strain change for dip angle in the range of 6-0°for A-2+ model with the same fault size, slip amount and strike angle as the original ones [fault geometry is listed in Table 5 of Ariyoshi et al. (2021)] by using a simple rectangular fault model (Okada, 1992) in a homogeneous elastic medium with a Poisson's ratio of 0.35 (Araki et al., 2017) and a rigidity of 20 GPa. The result shows that the volumetric strain at C0002 due to the fault model of A-1 is nearly constant (0.001 μstrain) and largely independent of its dip angle. Hereafter, we focus on A-2+ model because the dip angle lower than 6°for the shallower fault (A-2+ model) is more reasonable than that for the deeper fault (A-1 model).
Table 2 and Figure 10 shows the volumetric strain changes of Pc1 at three boreholes as a function of dip angle for A-2+ model in Figure 1A. If the dip angle of A-2+ is treated as 6°, the volumetric strain at C0006 is slightly negative, which is inconsistent with the sense of the observational result (+0.01 μstrain) in Figure 9. If we adopt the reasonable range of the volumetric strain change at C0006 (vertical axis) as between 0.006 and 0.014 in Figure 10 (vertical axis), this result suggests that possible dipping angle of A-2+ is lower than 4°, which is largely consistent with the volumetric strain changes at C0002 (−0.074 μstrain) and C0010 (−0.22 μstrain)  by applying 5.7 and 4.7 kPa/μstrain, respectively (Davis et al., 2009;Wallace et al., 2016;Araki et al., 2017). This result is just a simplified interpretation and the unified fault model is to be proposed in our future study.

Oceanic Fluctuation Impact on the Seafloor
To detect SSEs from pore pressure recorded in real time, it is important to distinguish oceanic perturbation from crustal deformation and understand the relationship between them. Gomberg et al. (2020) pointed out that ocean-pressure fluctuations on timescales longer than ∼5 days may affect only the termination and not the initiation of SSEs in the Hikurangi subduction zone. Since the relationship between termination time and ocean fluctuation has not been investigated yet at the Nankai Trough, we discuss it in the context of our results focusing on the SSE in March 2020.
To examine this suggestion, we calculate pressure fluctuations of the ocean variation on the seafloor at 1-h intervals by integrating the seawater density based on JCOPE ocean data assimilation system (JCOPE-T DA) (Miyazawa et al., 2021) from the sea surface to the seafloor and adding atmospheric sea level pressure. Note that tidal fluctuations on timescales shorter than ∼5 days are not accounted in the assimilation system. Another advantage of using JCOPE-T DA is the fact that the observed seafloor pressure data detided by the theoretical tide has still fluctuation about ±5 hPa and seems not reflect well baroclinic responses driven by Kuroshio meander . Figure 11 shows the time history of pressure fluctuated by ocean loading and pore pressures at Pc1 in each borehole by applying (Pc0, Pc0, KMC11) to the reference seafloor pressure at (C0002, C0010, C0006), respectively, with a time lag between the variations at C0006 and KMC11 listed in Table 1. To eliminate the tidal effect possibly included in the pressure fluctuation of the JCOPE-T DA, we adopt moving average for 25 h as shown by curves in Figure 11A. All the time histories of the moving average seem to represent short-term (about 5 days) oscillation. On the Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 717696 termination of the SSE, considering that fault regions of models A-1 and A-2+ cover under the boreholes at C0002 and C0010 (Figure 1), we find that their short-term oscillations are in increase phase. The amplitude of the short-term oscillation at C0002 and C0010 on the termination of SSE is about 6 hPa, which is consistent with statistically significance for pressure distributions at the SSEs in the Hikurangi as shown by Figure 3 of Gomberg et al. (2020). For C0006, Figure 11A shows the short-term oscillation is in increase phase with small amplitude (1∼3 hPa). Gomberg et al. (2020) suggested the physical relationship between the ocean loading and the SSE termination (in Figure 5 of their paper). The ocean loading affects the termination of the SSE as well as initiation when the fault slip causes temporal weakening of the fault strength. The fault strength is thought to depend on the time history of slip velocity and shear/ normal stress in addition to frictional properties (e.g., Yoshida et al., 2020) perturbed by the SSE. In order to estimate the time history of the fault strength, it is necessary to know the time history of detailed slip process in realistic model of the SSE as well as the frictional properties, which is our future work.
The seafloor pressure changes between the onset and termination of the SSE as indicated by long arrows in Figure 11A seem to be negatively correlated with the pore pressure change in Figure 11, particularly clear for C0002 (∼+6 hPa) and C0006 (∼−6 hPa). For C0002, the pore pressure change seems to be fitted as a ramp function, which is similar to the previous SSE in 2015 (Araki et al., 2017). By applying regression analysis, we estimate the time history of crustal deformation as shown by light cyan lines. For C0010, it is difficult to find the negative correlation because the crustal deformation component of the pore pressure is complicated due to the passage of fault slip by models A-1 to A-2+ in Figure 1A, which is statistically estimated by Ariyoshi et al. (2021) as shown by arrows in Figure 11B.
During the SSE in March 2020, a Kuroshio meander occurred around DONET stations. To trace the Kuroshio meander and to monitor the deep thermal structure of the cold (higher density) Target region is the same as Figure 1A.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 717696 core eddy south of Enshu-nada, we investigate the distributions in water temperatures at depths of 200 and 600 m (Figure 12), respectively. Associated with a westward movement of the current path of the Kuroshio indicated by the 15°C isotherm at a 200 m depth, a warm (lower density) water mass proceeded from the east of C0006 to near C0002 ( Figures 12A-C).
Meanwhile, at a 600 m depth, the deep cold water mass beneath the cold eddy off Enshu-nada propagated northward ( Figures 12D-F). This feature is typically known as S-shaped meandering pattern ( Figure 13 of Kasai et al., 1993). It is well known that the Kuroshio front variability includes 5-8, and 10-12 days periodic fluctuations , and 20-50 days fluctuations as S-shaped meander (Kasai et al., 1993;Miyama and Miyazawa 2014), where the transition between warm water intrusion and cold water detour seems S-shaped ( Figure 12G). The seafloor pressure changes (Eq. 4) are consistent with the variations of water density changes associated with those in temperature (the third term of the right-hand side of Eq. 4), and contributions from the other terms (the first and second terms of the right-hand of Eq. 4) are negligible. The density changes induced by the movement of the cold water are more evident at C0002 and C0010 than at C0006 (Figure 12).
Since less normal stress on the fault promotes slip generation on the basis of Coulomb Failure Stress (CFS) ( τ ＋ μ(σ-P), where τ is shear stress, μ is friction, σ is normal stress, P is pore pressure) (e.g., Labuz and Zang, 2012), we discuss the possibility of SSE triggered by ocean fluctuation quantitatively. According to Kanamori and Anderson (1975), stress drop is obtained as where Δτ is the stress drop, G is the rigidity, D is the slip amount on the fault, W is the fault width. For model A-2+, Δτ is estimated as 24 kPa, where G is treated as 20 GPa. Neglecting dip angle (<4°), we can roughly estimate ΔCFS [ μ(Δσ), where we neglect dip angle <4°] as 2∼3 hPa, which is about 1% of the stress drop. From theoretical research, static stress change is more effective to trigger fault slip than dynamic stress perturbation (e.g., Gomberg et al., 1998). Yoshida et al. (2020) performed a numerical simulation of earthquake cycle on the basis of rateand state-friction laws on a simplified subduction plate boundary. With static and dynamic change of effective normal stress 1/8 time of recurrence interval before the earthquake onset, earthquakes can be instantaneously triggered by static change of less than 2.5% of stress drop while not by dynamic change of greater than 8% of stress drop. If we assume that the frictional state on the fault model A-2+ as the condition ready to occur fault slip due to nearby slip on the fault model A-1, and that the trend change is approximately applicable to static change, the amount of pressure for the trend change could be sufficient to promote slip upward through the source regions of sVLFEs.
Focusing on the pressure change around the termination of the SSE in Figure 11B (Supplementary Figure S1), impulsive pore pressure change occurs all the three borehole stations. This is thought to be excited by seismic waves of the March 25, 2020 Mw7.5 earthquake in northern Kuril Island (Ye et al., 2021). The arrival time of the P-wave is estimated as 298 s later after the origin time (11:49:21 in JST + 9 h from UTC) from tauP (Crotwell et al., 1999), which seems consistent with the arrival time at KMD13 (11:54:19) and later than the SSE termination time at around 08:40 in JST (Supplementary Figure S1).
Since the standard deviation of the estimated SSE onset/ termination time is 9 h , it might be possible that the seismic waves also affect the termination. Since the previous studies focused on the clock advance of triggering earthquakes (Gomberg et al., 1998;Yoshida et al., 2020), it is necessary to investigate the clock advance of terminating earthquakes, which is our future study. Figure 13 shows a schematic explanation about the possible scenario of the SSE in March 2020. The trend change driven by oceanic fluctuation (Figures 13A,B) promotes shallower slip causing sVLFEs during the passage of Kuroshio meander, and FIGURE 13 | Schematic scenario of the possible process for the SSE with ocean loading (A) during model A-1, (B) during model A-2+, and (C) termination of the SSE. Representative focal mechanism of VLFE is obtained from Ariyoshi et al. (2021). Double arrows represent volumetric strain change due to crustal deformation driven by the SSE.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 717696 increase of normal stress due to the short-term oscillation ( Figure 13C) promotes termination of the SSE as observed in the Hikurangi subduction. Therefore, monitoring of oceanic oscillation especially ocean current meander is also important to forecast the growth of SSE and avoid the misdetection of crustal deformation. Because JCOPE-T DA is the latest version of JCOPE (Miyazawa et al., 2021), it costs much time to get the reanalyzed ocean modeling going back from the present to past time. It is our future work to apply JCOPE-T DA to the time when previous SSEs occurred, which would reveal more detailed relationship between Kuroshio meander and previous SSEs.

CONCLUSION
Our conclusions are summarized as follows.
i) If we use the reference pressure on the seafloor just above the borehole, time lag is negligible for one second sampling. This shows we can monitor pore pressure simply by subtracting the reference pressure with valid factor if the instrumental drift component is negligible. ii) We demonstrate that the reference seafloor pressures at boreholes can be substituted with nearby DONET seafloor pressure gauge data. This would give more robust condition to monitor the pore pressure in case of malfunctioning of the borehole reference pressure gauge. iii) The time lag of the pressure between borehole and nearby DONET station is temporarily perturbed due to tidal oscillation, which can be explained by superposition of tidal modes. The amplification factor of the reference seafloor pressure to borehole pressure is perturbed during the occurrence of nearby SSE in case of significant pore pressure change. iv) Using the seafloor pressures recorded at the nearby DONET station to the borehole pressures at C0006, we can extract the pore pressure change due to crustal deformation of the SSE in March 2020 at the deepest borehole (Pc1) of C0006 about 0.6 hPa. v) On the basis of ocean modeling, the ocean impact of seafloor pressure without tidal effect has short-term (5 days) oscillation and trend change between the onset and termination of the SSE. The amount of ΔCFS around C0006 is about several percent of stress drop for the SSE, which may be sufficient to promote updip slip propagation toward the fault model which triggered sVLFEs.
In summary, real-time monitoring of precise pore pressure with nearby DONET stations and ocean modeling would be a powerful tool to judge the onset and termination of SSEs promptly.

DATA AVAILABILITY STATEMENT
Borehole data is available from the website of the J-SEIS (JAMSTEC Ocean-bottom Seismology Database; https://join-web.jamstec.go.jp/ join-portal/). DONET data is available from the website of NIED (National Research Institute for Earth Science and Disaster Resilience) https://doi.org/10.17598/nied.0008. The JCOPE-T DA data are available from the website of JAMSTEC; http://www. jamstec.go.jp/jcope/htdocs/e/distribution/.

AUTHOR CONTRIBUTIONS
KA estimated the best solution of the factors and time delays between boreholes and nearby DONET station, in addition to write the manuscript overall. TK tested the application of nearby DONET seafloor pressure data. YM wrote Ocean model product used for calculation of atmosphere-ocean-induced sea floor pressure fluctuations section. SV customized the JCOPE-T DA for boreholes and DONET stations. TI calculate the dependency of the volumetric strain on the dip angle in Figure 10. AN revised Oceanic fluctuation impact on the seafloor section. JG revised the manuscript as a native English writer. EA designed the borehole and DONET observatories. TM advised the interpretation of the results from JCOPE-T DA. KS managed the borehole and DONET data server. SY supported for making Figure 13. TH and SK coordinated the research discussion. NT managed the data sharing of DONET.

FUNDING
This work was partly supported by the Japan Society for the Promotion of Sciences (JSPS), Grant-in-Aid for Scientific Research (KAKENHI) Grant Numbers JP20H02236 and JP20KK0097.