Overpressure Estimation Based on Velocity Data Analysis in The Subduction Zone of Hyuga-nada at Western Nankai Trough, Japan

Subduction zone is one of the areas which has great potential of earthquake. Nankai trough in Japan is such an area with subduction zone where earthquake records are well archived, despite its complexity. Also, slow earthquakes have been detected in this zone over the past 20 years, so that there are many opportunities to understand the dynamics and deformation in such an active setting. To expand the knowledge about slow earthquake which might occur and trigger a bigger earthquake, this study wants to infer the relationship between overpressure, one of the states and properties of seismogenic zone, and tremors. To carry out this research, 4 velocity models of seismicity line taken in Hyuga-nada are used. And the recorded tremor data placed around the lines came from Yamashita. In correlating the overpressure and tremors, this research is started by converting the Vp model to porosity model using Wyllie equation to predict the pore pressure conditions in the area. Then, the model will be used for fitting the actual data. Assuming there is an exact location following the normal compaction curve as a function of effective stress, it can be used for estimating the overpressure. Overpressure is economized by comparing the observed porosity and the reference curve at the same depth below seafloor. Potential overpressure occurred may be associated with areas of low Vp. The estimated overpressure is in the range of 1.2 to 12 MPa in the depth range of 700 to 4000mbsf. This study also concludes that the overpressure zone overlaps with the tremor location, although the overpressure value estimated is also still uncertain due to the uncertainty of the density value used.


Introduction
Earthquakes can cause liquefaction, landslides, fires, and tsunamis that will lead to significant levels of damage and loss.Most earthquakes occur at boundaries where plates converge, diverge, or move laterally past each other.Convergent boundaries can also occur in regions of continental collision that produce tectonic compression [1].Regions with complex characteristics such as subduction zones provide an opportunity to study the interaction of deformation and fluid flow in active settings.Fluids in subduction complexes play an important role in deformation, geochemical and thermal transport, and seismogenic zone behavior [2].However, it is not well known whether pore fluid pressure is localized within the fault zone or permeates the entire forearc area [3].While it is known that pore fluid pressure plays a key role in the generation of subduction zone earthquakes, it is unfortunate that the quantification 1307 (2024) 012023 IOP Publishing doi:10.1088/1755-1315/1307/1/012023 2 of pore fluid pressure around subduction zone faults has, to date, been limited [4].Only 2% over the last 20 years.In fact, it is not only megathrust earthquakes that occur in complex subduction zones, but there are known to be slower earthquakes that occur in areas near the boundaries of updip and downdip zones in locked interseismic and cooseismic slip.These earthquakes are said to be associated with excess pore fluid pressure that triggers a decrease in effective stress [5].One of the most complex subduction zones is the Nankai accretionary complex; formed by the subduction of the Philippine Sea Plate beneath the Eurasian Plate in southwestern Japan which is creating a convergent plate boundary.There are also several seamounts rising from the seafloor in this area.The subduction zone in the Nankai Trough is a typical subduction system characterized by the subduction of several geological units of the Philippine Sea Plate or PSP.These units include the Kyushu-Palau Ridge, Shikoku Basin, Kinan Seamount Chain, and Izu-Bonin Arc.The PSP is a rhombus-shaped plate with a narrow NS direction, which is surrounded by several convergent plate boundaries [6].Sediments scraped from the top of the PSP or the subducted Philippine Sea Plate have formed a 5 -10 km thick accretionary prism landward of the Nankai Trough (Moore et al., 1990).In 1944, the 100-km-wide locked zone of the subducting plate interface shifted, causing the M8 Tonankai earthquake that was destructive and even devastating [7].Very low frequency (VLF) earthquakes have also been approximately detected near the trench axis of the Nankai Trough with a dominant frequency of 0.1 Hz.It has also been discussed that VLF earthquakes may occur as a result of slip on the fault plane in the mega-splay fault that branches off from the seismically coupled zone of the megathrust earthquake [5].Meanwhile, Obara (2016) mentioned that not only VLF earthquakes appear based on the Nankai subduction cross section.There are various types of slow earthquakes except VLF; slow slip events (SSE), episodic tremor and slip (ETS) [8], low frequency earthquake (LFE) [9].To investigate the occurrence of slow slip earthquakes; events in distribution zones associated with strike heterogeneity at different depth ranges, this study intends to examine the Hyuga-nada area in the western part of the Nankai Trough to perform pore pressure prediction to answer the relationship between slow earthquakes and Low-Vp areas.The author designed this research because it is known that there are many previous studies that have studied and provided insights into the state and properties of the Nankai Trough seismogenic zone [4][7] [10].One of them has detected so-called slow earthquakes, slow slip earthquakes, in this subduction zone over the past 20 years.Type of slow earthquakes such as very low frequency earthquake and tremor.The authors hope to provide conclusions from the relationship between overpressure and tremor recorded there.

Methodology
To address the above issues, this study uses velocity data modeled from extracted seismic reflection data taken at Hyuga-nada, western Nankai Trough.A total of 4 seismic trajectories were extracted into 2dimensional compressional velocity model for each line.Then, as a correlation analysis material in this study, the slow earthquake data used is the tremor type taken from Yamashita (2021) for 5 years of tremor data from 2013 to 2018 around Hyuga-nada.  2.1.Converting Vp to Porosity Screaton (2002) in his study, examined the evolution of porosity near the Nankai accretionary prism and looked at the implications for pore fluid pressure.Porosity was found to decrease with depth, and to have a porosity-compactness profile in accordance with the definition from [15].Therefore, this study first converts the compressional velocity or Vp value into a 2-dimensional model of porosity.Through nonlinear bi-square curve fitting, the Wyllie (1959) equation was calculated specifically for the Nankai Trough region as follows [16]: The Wyllie equation was chosen because it is a global equation that is still commonly used to cross-plot velocity data with derived porosity and resistivity.It has the advantage that it can also be forwarded for conversion into 3D porosity volumes.In addition, compared to other formulas such as Hoffman's, Hyndman's, Erickson & Jarrad's, the porosity profile calculated using Wyllie seems to be more in line with the conditions in the Nankai Trough.Therefore, this function was applied to convert the porosity values from the Vp model data of each line, namely HYU01, HYU02, HYU18, and HYU22, so as to obtain 4 2-dimensional porosity profiles with respect to depth.

Effective Stress Calculation
As shown in Figure 5, one input data Vp will get two treatments: converted to porosity by Wyllie and created theoretical curve of porosity.In this case, the porosity depends on the effective stress.Meanwhile, the effective stress is equal to the result of reducing all charges above the depth of the area of interests by hydrostatic pressure and pore pressure.The relationship is as mathematically written as follows: ) where   is pore pressure,  ℎ is hydrostatic pressure, and ∆P is the overpressure to be seen.
The effective stress is also correlated with the physical property of density, where the function that states the effective stress can be calculated by integrating the density value with depth is as follows: ( ′ ) =  ∫   ( ′ ) ′ +      0 (4)  ℎ () =  ∫   ()  0 (5) where   ( ′ )is the density of the bulk,   is the density of seawater, and  ′ is the depth of the water.

Making a reference porosity curve
Because of the limitation of drill log data, in the calculation of hydrostatic pressure, the density value obtained from the reference sediment density for   ( ′ ).And for the   , we used the value of the static density of water.Here, the hydrostatic pressure represents the pore pressure value, because in this theoretical porosity model is considered a normal condition with no overpressure or we can state that the ∆ = 0.This is the nature of porosity which the properties assumed to correspond to its character; decrease exponentially as the equation defined by Rubey  (10) The reference porosity profile with respect to depth is obtained after interpolation and finding the value of the  0 and .This step is repeated for a distance of different horizontal columns.These differences are then considered in order to select the best reference point.Since this process is done on each line, several values of  0 and  will be obtained for each line.
To select the best reference in this stage, we calculate the mean absolute error (MAE) for each selected column point.This error calculation is carried out with the concept of finding the difference in values from two sets of data that are observed to have similar conditions or phenomena that applied [11].The observed data sets here are the porosity data converted using Wyllie, and the reference porosity data and the porosity reference data calculated using the stress theoretical formula.Mathematically, the process can be written as follows: (11) where ei is the absolute error, yi is the calculation data, and xi is the actual known data.In this context, yi is the reference porosity and xi is the porosity based on Wyllie.The horizontal x position which is chosen as the reference will affect the estimation value of overpressure calculation.The difference in the reference horizontal x will play a role in shifting the position or anomalous layout trend.In addition, it will also affect the size of the estimated value due to the slope variables (), , and  0 that comes from this step.

Estimating Overpressure
On the next step, the overpressure will be estimated.So, different with the best-fit calculation where the condition is assumed as normal with no overpressure, the calculation in this step will consider the overpressure is in there.The following equation is used to estimate the overpressure value:  =  0  −  ′ ,   ′ =   − ( ℎ + ∆) Overpressure is obtained by comparing conditions at the same depth for two data sets: those that are assumed to be normal (reference) and those that are as they are.For each column, the distance x horizontal is estimated, and a map of its distribution is made in 2D with depth.With complete variable value needed, this process is done using MATLAB.
Finally, to analyze the correlation between the results of the overpressure estimation and the occurrence of tremor in the study area, the pore pressure points with overpressure conditions were plotted to the area on the siesmic trajectory using QGIS, and their positions were reviewed.To display the final results, the anomalies were mapped on the bathymetric map with the tremor and the seismic lines in a figure to analyze it and find the correlation.

Result and Discussion
The initial input data, Vp 2D, has vary value with depth.After converting the Vp to porosity, the value trends have a similar shape, exponential to depth, following the nature of its characteristics.Both the Vp and porosity are inversely proportional to depth.The conversion result as shown in Figure 7b with initial Vp model as shown in Figure 7a.
Visualization of the distribution of variations in Vp and porosity data as shown in the figure indicates the depth of seafloor from the sealevel.This display will give an idea of zones that are normal without anomalies or not with overpressure anomalies.The pattern of value distribution and trajectory direction were considered by the author in choosing the reference point or data column before making the reference curve.Sampling of several horizontal x points is performed to find the best model that most closely approximates the best-fit of its curve.
The horizontal x chosen as a reference is best, represented by its small error value.The difference in choosing x horizontal references will play a role in shifting the position or anomalous layout trend.In addition, it will also affect the size of the estimated value.It is because of the slope, compressibility, and porosity variables will be obtained zero in this stage.Then, to obtain these values, density is needed in the calculation using a formula derived from Athy (1930) as written in Equation ( 7).This study took a reference density value from [12] which is 1700 kg/m3 for solid material, and for static water we put 1300 kg/m3.The result of determining each x reference as shown as black line in figure 7b, is presented on Table 1 including the variable needed and the minimum error.The results in the table will then be used as a basis for calculating the overpressure at the depth of each line.The overpressure value is obtained by comparing conditions at the same depth for the two data sets.The estimated value is gradually calculated per column of porosity data and presents the results as in Figures 8-11.The interpolation corresponds to the x position at the distance of each selected trajectory, as indicated by the black line in Figure 7b (red line in the overpressure distribution map).Figure 8 shows the presence of overpressure conditions at the depth of HYU01 traverse as marked by the white box.In the figure, the seabed began to be identified at depths of 2000-3000m below sea level.In the HYU01 trajectory, anomalies are seen at depths of 4500m and 6000m from sea level or 2500m below seafloor (bsf) and 4000mbsf with anomalous values of 4 MPa and 1.2 to 3 MPa, respectively.Estimation based on reference point at x=60550m with error 0.028.On the beginning of seafloor, hereinafter referred to as the shallow area, the porosity value ranges from 0.55 decreasing to 0.35 at depth.Then, in the medium area to a depth of 10000m, the porosity value is around 0.3.And at more depths than that it shows a decrease in the porosity value again to less than 0.1.As can be seen again in Figure 7a, this plot evidently depicts the area of higher velocity anomaly which correlates to a lower porosity zone.
In Figure 9, the seafloor in HYU02 appears deeper than HYU01 and the anomaly is shown at a depth of 5000-7000m from the surface or 2000-2600mbsf with an indication of a larger overpressure value, estimated at 6-12 MPa.This value was obtained by taking a reference site at x=45000 with an error of 0.024.The trajectory of HYU02 is visualized in Figure 9, with the seabed jutting deeper and deeper from a depth of 2400m to 5000m below sea level.The velocity anomaly stands out, as seen from the velocity immediately increasing at a depth of 6000m at a horizontal distance of 20000-40000m trajectory.It is referred to as an anomaly because at horizontal distances of 0-10000m and 41000onwards it has not experienced an increase in velocity at the same depth.And the porosity value profile follows, making it look like a mountain of mountains is poking down from below.Porosity is very small at a small horizontal distance since in areas that tend to be shallow, which is 0.1.But the important factor is the small area of lower velocity below the higher velocity anomaly.Lower velocity correlates to a region of higher porosity which suggests that the pore fluids support a disproportionately large part of the leading of the overpressure.
Showing in Figure 10, the HYU18's seafloor began to be identified at a depth of 2500m below sea level.This figure also shows that in the seafloor area, the porosity value ranges from 0.6 decreasing to 0.5 at depth.Then, in the medium area to a depth of 10000m, the porosity value is around the figure 0.2-0.3.Deeper than 10000m, the porosity value continues to decrease and tends to be constant from a depth of 15000m below sea level.There is a shallow anomaly with an overpressure value of 5-6 MPa at a depth of 3200m from the surface or an estimated 700mbsf.The reference site for this HYU18 traverse is x=9300m, with an absolute error of 0.024.
The HYU22 trajectory is plotted in Figure 11, where the seafloor begins to appear at a depth of 2000m below sea level.With the trend of forming two mountains, the velocity continues to increase until it tends to be constant.The improvement is regular and tight.Once converted, it can be seen in the shallow area of the seafloor, the porosity value ranges from 0.6 decreasing to 0.5 on the seafloor face along the track.Then, in medium to 10000m deep areas the porosity value ranges from 0.25-0.3,but at a horizontal distance of 60000m the value has decreased to 0.1.To a depth of more than 15000m, the porosity value decreases until it tends to be constant.Next, the figure shows the overpressure distribution profile with the anomalous trend by calculating using the selected reference site at x=54300m on the HYU22 line.The overpressure value obtained is marked with a white box, ranging from 1MPa at a depth of 5000m from the water table or equivalent to 2500mbsf.Note that this result was obtained with an average error of 0.017.Overpressure zones indicate restricted vertical and horizontal permeability, possibly by a high density of faults growth, high shale-to-sand ratios, deposition of impermeable sediments, and montmorillonite diagenesis among other factors [13].Overpressure might be there, is caused by fluid pressure in the rock body or it might also have correlation with temperature at depth which can lead to dehydration that helped the underthrust sediments to retain such high porosities.This might be a higher chance of slips and or earthquake can occur.So, to see the distribution of overpressure above the previously recorded tremor event in the Hyuga-nada area, the marked anomalies were plotted on the surface of the bathymetry map that had been overlaid with the data tremor coordinates (see figure 4).The result is as shown in Figure 12.
Tremors recorded a total of 6888 events.For each line, approximately crossed 872 tremors at a distance around HYU01, 287 tremors in the HYU02 line's area, 156 tremor events passed through the HYU18 line, and as many as 421 tremors crossed by the HYU22 line.Four histograms are made to see whether the overpressure estimated has correlation with the tremor peak area or not.To simplify the visualization, Figure 12 displays the orange bullets overlayed with the peak of histogram regarding the thickness of the black dots which represents the tremor frequency around Hyuga-nada. Porosity anomaly can be converted to overpressure, which is both a function of space and time.It is usually still difficult and challenging, partly because of the complex geology and the porosity trend varies in place.

Conclusions
The estimation of porosity, pore pressure, and overpressure values in Hyuga-nada western Nankai trough, Japan, was carried out by applying Wylie Equation to the Vp mode, followed by determining the reference location assuming that the normal compaction zone is exponential as a function of effective stress.The high pore pressure or known as overpressure was estimated by comparing the observed porosity with the reference curve at the same depth.From the four lines studied, overpressure potentially exists in the low Vp area, with an estimated magnitude of 1.2 to 12 MPa in the depth range of 800mbsf to 4000 mbsf.Then, the estimated overpressure zone can be inferred to overlap with the points of occurrence of tremors scattered in Hyuga-nada.This is a possible correlation between the overpressure potential and the distribution of tremors in western Nankai Trough, Japan.
Meanwhile, there are still doubts about the accuracy of the choice of the reference sites, whether the method for choosing the reference is the best and can be applied to all conditions or not.Therefore, a better method or technology should be used to determine the reference site.In addition, the estimated overpressure value also needs to be reconfirmed because the calculation is estimated with a density reference value not data value taken directly from the research field.It is necessary to support density log data to make this research mor advance.

Figure 1 .
Figure 1.Illustration of the spatial variations in the frictional properties of a subduction fault: (a) The schematic characterization of the megathrust frictional environment, showing 'seismic' region's as frictional sliding and regions labeled as 'aseismic' stable or episodic slow.'Conditionally stable' labeled area which displace aseismically except when accelerated by failure of adjacent seismic patches.(b) The aseismic updip zone overlied by the outerwedge, and the seismogenic zone overlied by the inner wedge.The shallowest part of the subduction has its strength of interface extending to several depths increases during an earthquake, resisting rupture [14][18].

Figure 3 .
Figure 3. Types of slow earthquakes at the plate interface of the Nankai Trough subduction zone [8]

Figure 4 .
Figure 4. Tremor distribution and seismic lineaments at Hyuga-nada used.The black dots are the tremors, and the yellow lines are the seismic lines.

Figure 6 .
Figure 6.Three random x column in the same seismic line, display the difference between actual data and its best-fit data porosity profiles in Ms. Excel.The orange profiles the actual data based on Wyllie, and the blue line profiles the best-fit data for the same x column.The best for the reference might be the most similar one, it will have the minimum error

Figure 7 .
The visualization for 2D Vp model in m/s (a), the converted porosity map using Wyllie, with the black line is the x reference chosen (b), and the 1D curve of the x reference data at depth (blue curve) with its best-fit curve based on calculation (red curve) (c) IOP Publishing doi:10.1088/1755-1315/1307/1/0120238

Figure 8 .Figure 9 .Figure 10 .Figure 11 .
Figure 8. Porosity map (left) and the overpressure estimation for line HYU01 with the reference line (black-red) at x=60550m show the anomaly in around 3MPa in the white box

Figure 12 .
Figure 12.Overlayed map of overpressure and tremor distribution from the surface.The oval marks the anomaly position from each line processing

Table 1 .
Position of x reference for each line and the obtained value of ,  0 , , and error calculation