Identification of active faults and tectonic features through heat flow distribution in the Nankai Trough, Japan, based on high-resolution velocity-estimated bottom-simulating reflector depths

Estimates of heat flow can contribute to our understanding of geological structures in plate convergent zones that produce great earthquakes. We applied automated velocity analysis to obtain the accurate seismic profiles needed for precise heat flow estimates using six new seismic profiles acquired during R/V Kaimei KM18-10 voyage in 2018. We calculated heat flow values in the accretionary wedge of the Nankai Trough off the Kii Peninsula, Japan, from the positions of widespread bottom-simulating reflectors (BSRs) in seismic reflection profiles. Calculated conductive heat flow values from the depth of the BSR agree with previous studies where a regional trend is observed from ~ 50 mW/m2 to < 40 mW/m2 60 km landward from the deformation front. This trend is caused by thickening of accretionary sediments and the subduction of the Philippines Sea plate. Segments of profiles are marked by anomalous high heat flow values. Such anomalies represent alterations of the shallow crustal thermal structure caused either by a combination of topographic affects, surface erosion of the seafloor, or by fluid flow that transports heat by advection. We interpret heat flow anomalies (~ 100 mW/m2) as indicators of active faulting, which correspond to low seismic velocity zones along faults. Our results also showed relatively high heat flow at the landward end of several survey lines close to the Kii Peninsula, which we interpret to the possible presence of plutonic rocks that underlie the Kii Peninsula and extend offshore and may be the cause of geothermal springs, steep geothermal gradients, and high heat flow.


Introduction
Heat flow data are useful for interpreting geological features and process in subduction zones.Since heat flow can be affected by several different processes (e.g., the ascent of fluids toward the seafloor, rapid erosion and sedimentation, and tectonic uplift), we can infer aspects of tectonic activity or thermal regimes in subduction zones from heat flow data.Heat flow distribution has been investigated by several methods, including the use of seafloor sediment probes (Christoffel and Calhaem 1969), piston core samplers (Kinoshita et al. 2008), and the depth distribution of methane hydrate deposits derived from the position of bottom-simulating reflectors (BSRs) in seismic profiles (e.g., Yamano et al. 1982;Davis et al. 1990;Hyndman and Spence 1992).
BSRs, which mark the phase boundary between free methane gas below and methane hydrates above (Markl et al. 1970), appear in seismic profiles as high-amplitude reflectors with reverse polarity compared to the seafloor reflections and are nearly parallel to the seafloor (Vanneste et al. 2001).BSRs can be disrupted by several processes that also affect heat flow, allowing us to infer thermal conditions in the shallow sediments where BSRs appear.
In this study, we derived estimates of heat flow from the distribution of BSRs in the Nankai Trough and used this information to investigate geological features in the accretionary wedge.This method can provide continuously distributed heat flow data along seismic transects with finer resolution than direct in-situ measurement methods, which may help reveal geological structures.However, estimating heat flow from BSR profiles requires accurate depth information.In this study, we applied automated velocity analysis (Fomel et al. 2013) to obtain a high-resolution P-wave velocity structure (Mukumoto et al. 2019;Kret et al. 2020), which helped reveal the presence of free gas beneath the methane hydrate layer as low-velocity anomalies (Chhun et al. 2018;Eng and Tsuji 2019).Automated velocity analysis can incorporate all available seismic data to produce a seismic velocity model based on dense spatial coverage and hence yield seismic profiles with improved accuracy.In this paper, we describe our estimates of conductive heat flow derived from the improved depth model and discuss the geological structures inferred from the estimated heat flow anomalies.
The Nankai Trough is a well-studied plate subduction zone that has been investigated with a variety of

Graphical abstract
geophysical techniques including seismic, gravity, and electrical resistivity surveys (e.g., Honda and Kono 2005;Kasaya et al. 2005;Kodaira et al. 2007).Seismic reflection profiling is one of the most effective of these for characterizing the details of geological structures and physical properties.The detailed knowledge of heat flow derived from the high-resolution velocity and reflection profiles can help identify active faults that may affect thermal conditions.These inferred active faults could inform further geophysical surveys, for purposes such as optimal deployment of seismometer networks to monitor fluid flow and fault dynamics.

Geological setting
The Nankai Trough is a convergent plate boundary where the Philippine Sea plate is subducting beneath the Eurasian plate at approximately 3-5 cm/year (e.g., Seno et al. 1993).The Nankai Trough has been divided into several segments on the basis of the rupture zones of historical great earthquakes (e.g., Ando 1975;Cummins et al. 2002;Hori 2006), and one such segment is located off the Kii Peninsula.On the landward side of the trench, a sedimentary accretionary prism has formed on the Eurasian plate, overlain by a forearc basin.The imbricate thrust zone (ITZ), in the middle of the accretionary prism, is defined as weak folding and inclined turbidite beds consisting of fractured mudstone and deformed sand layers (Kawamura et al. 2009).
The accretionary prism accumulates methane driven by upward fluid migration, and in a steady-state condition, the BSRs correspond to the base of the gas hydrate stability field (BGHS) (Grevemeyer and Villinger 2001;Yamada et al. 2014), and these are widely distributed in the Nankai subduction zone (Baba and Yamada 2004;Colwell et al. 2004).The forearc basin contains methane hydrates within alternating sediment layers of sand and mud (Suzuki et al. 2009).Fluid containing methane that originates in the décollement migrates upward through faults and intergranular porosity during the accretion process (Yamada et al. 2014) and is evident as cold seeps in the seaward side of the outer wedge (Ashi 1996;Ashi et al. 2002a;Tsuji et al. 2014).
Heat flow has been extensively measured and estimated by BSRs offshore the eastern Kii Peninsula (e.g., Kinoshita et al. 2008;Hamamoto et al. 2011).However, heat flow data are sparse off the western part of the Kii Peninsula.Moreover, several geophysical studies, including seismic and gravity surveys, have suggested the presence of a pluton beneath the Kii Peninsula and its south offshore area (Kodaira et al. 2006;Kimura et al. 2014Kimura et al. , 2022;;Arnulf et al. 2022).Those studies have suggested that this plutonic rock intruded during the middle Miocene into the upper plate of the island arc crust.However, the relationship between this intrusion and the thermal structure of the Kii Peninsula is still unclear.

Seismic data
The Japan Agency for Marine-Earth Science and Technology (JAMSTEC) acquired two-dimensional multichannel seismic data and high-resolution bathymetric data collected from the central Nankai Trough in 2018 with R/V Kaimei (Nakamura et al. 2022).In this study, we used six survey lines from that survey to investigate the relationship between heat flow distribution and active faults (Fig. 1).A tuned airgun array of totaling 173 L with 13.8 MPa (2000 psi) was fired with a shot interval of 50 m, and reflected waveforms were recorded on a towed hydrophone streamer 4.5 km long with a channel interval of 3.125 m.Data were acquired at a sample interval of 2 ms and a record length of 14 s.

Seismic processing and automated velocity analysis
Before performing automated velocity analysis, several signal enhancement and noise reduction techniques were undertaken (Fig. 2).In this study, we used multi-channel seismic data after preprocessing, including noise reduction, deghosting, designature, debubble, and multiple reflection attenuation (Nakamura et al. 2022).In addition, we applied a bandpass filter with a frequency range of 5-90 Hz to extract the appropriate frequency range of dominant reflection signals.A predictive deconvolution filter was then applied to remove the effects of swell or bubbling.To reduce the computational cost, we decimated shot domain gathers from a 3.125 m to a 6.25 m interval, thus changing the common mid-point (CMP) interval to 3.125 m.After the sorting process, the total number of CMPs was approximately 24,000 for each line.
A high-resolution P-wave velocity structure, determined by velocity analysis, is crucial to accurately map BSR depth and estimate heat flow distributions.We applied automated velocity analysis for this purpose (Fomel et al. 2013), allowing us to estimate seismic velocity from all CMP gathers.Automated velocity picking takes place on calculated semblance panels (Fomel 2009).We first generate a semblance panel of root-mean-square velocity, ranging from 1400 m/s to 2000 m/s with an interval of 20 m/s at every time sample.After estimating root-mean-square velocity by automatically picking semblance peaks (Fomel 2009), the interval velocity was computed based on Dix's equation (Dix 1955).We achieved a high spatial resolution in the resultant velocity structure by estimating seismic velocity from all CMP gathers.Finally, post-stack time migration (post-STM) and depth conversion was applied to all profiles.

BSR picking and heat flow calculation
After the depth conversion of seismic profiles, we picked BSRs manually at every 50 CMP number, which is equivalent to approximately 156 m.The lateral resolution and resolvable thickness of layers on seismic profiles is related to the Fresnel Zone and law of quarter wavelength (e.g., Lindsey 1989;Sheriff 2002).This gives the first Fresnel Zone diameter of 200 m and minimum resolvable thickness of 10 m for the KM18-10 processed seismic data.The picked horizontal interval, which is shorter than the lateral resolution, is reasonable since migration reduces the Fresnel Zone radius.
Since a BSR coincides with the BGHS (Grevemeyer and Villinger 2001) in a static system, the required parameters for heat flow calculations can be inferred from BSR depth.We calculated heat flow based on Fourier's law: where Q is heat flow [mW/m 2 ], k is the bulk thermal conductivity of the sediment between the seafloor and BSR [W/mK], T is temperature [K], and Z is depth [km].The thermal conductivity can be calculated from P-wave velocity by using an empirical relationship (Duffaut et al. 2018); however, scattering of P-wave velocities may lead to significant errors in thermal conductivity.To remove the effects of scattering, we assumed a thermal conductivity of 1.2 W/mK by averaging the estimated velocity values between the seafloor and the BSRs.Furthermore, we estimated the geothermal gradient from BSR depth by assuming hydrostatic pressure conditions (Ashi and Taira 1993).
(1) We used the following method to estimate thermal gradients.The hydrostatic pressure at the BSR was estimated from the BSR depth on the seismic profile.We estimated the temperature of the BSR by reference to the methane hydrate phase diagram determined by Hyndman et al. (1992) from laboratory experiments, considering the pure water and pure methane.The seafloor temperature was reported as 2 °C by Integrated Ocean Drilling Program Expeditions 315 and 316 (Harris et al. 2011), and we applied this value to the parameter T seafloor .We derived the thermal gradient from the temperature and depth information.Heat flow at each CMP was then estimated by calculating the product of thermal gradient and thermal conductivity (Eq.1).

P-wave velocity structure
The P-wave velocity structure obtained from our automated velocity analysis revealed detailed geological features in our study area.Since seafloor multiple reflections obscure the deeper part of this structure, we only display P-wave velocity profiles shallower than 1 s in two-way travel time below the seafloor (Fig. 3).This was adequate to identify the BSR, which is located at ~ 0.6 s below the seafloor.By comparing the velocity model and the seismic profile, we determined that the interval velocity range in the sediments containing BSRs was approximately 1800-2300 m/s.The noteworthy feature of velocity analysis is the low-velocity zones (< 1600 m/s), apparent on each velocity profile (Fig. 3, red arrows).The low-velocity zones coincide with major faults (e.g., out-of-sequence thrusts) coring anticline ridges (Fig. 3d) and also observed beneath BSRs (Fig. 3c).Low velocity zones are inferred to reflect the presence of free gas as reported by Chhun et al. (2018) in the Kumano basin and gas charged fluids migrating along faults into the hydrate stability zone (e.g., Korenaga et al. 1997;Tinivella and Accaino 2000).

BSR picking
We can identify the wide distribution of BSRs clearly on the migrated seismic profiles (Fig. 4).Spatial discontinuities of the BSRs are found around the thrust faulting in the accretionary prisms (Fig. 5).BSRs directly under the seafloor where sedimentation or erosion has occurred recently are not parallel to the seafloor or have an edge that diverges from parallel (Fig. 5).Indeed, we observed the double BSR on the seaward part of KIN_054 (Fig. 5b).We picked the lower BSR by considering the erosion effect as we discussed later.

Heat flow
Figure 4 also shows the estimated heat flow and the corresponding seismic depth profiles.The heat flow estimated from the upper part of the accretionary prism to the ITZ ranges from 20 to 100 mW/m 2 .From our results, heat flow typically increases in the trenchward direction (Fig. 6b), except at several points inferred as extremely Fig. 2 The flowchart of seismic data processing conducted in this study.Automated velocity analysis was applied on CMP sorted gathers after preprocessing (Nakamura et al. 2022) and additional signal enhancement.Final reflection profile was obtained through depth conversion of the post-stack time migration based on the high-density velocity model high anomalies (e.g., Fig. 4c, d, f ).Each seismic profile indicates multiple lateral trends in heat flow as well as extremely high heat flow anomalies as great as 90 mW/ m 2 (Figs.4d, 5a).
The maps of estimated BSR temperatures and heat flow values show a trenchward increase on all six survey lines (Fig. 6a, b).The range of BSR temperatures in our study area is between 17 °C and 25 °C (Fig. 6a), which is consistent with in-situ temperature records determined by long-term borehole monitoring in the Kumano forearc basin (Sugihara et al. 2014).The map in Fig. 6 is limited to the range of 30-60 mW/m 2 to emphasize local heat flow trends (see Additional file 1: Figure S1 for the map showing the full-width color scale of heat flow).Our heat flow distribution (Fig. 6b) generally agrees with previous studies (e.g., Marcaillou et al. 2012;Ohde et al. 2018) except near the eastern end, where our estimate was slightly higher (line KIN_070 in Fig. 6b).
Here, we note the three main sources of uncertainty in our heat flow estimates: (1) seismic velocity, (2) thermal conductivity, and (3) depth error derived from migration processing approaches.Although there are other possible sources of heat flow uncertainties, such as gas composition, salinity, seafloor temperature, and sedimentation, we did not include these factors because we could not confidently quantify their uncertainty.We hence evaluate several possible sources of uncertainty that could be estimated quantitatively: seismic velocity, assumed thermal conductivity, and picked BSR depth.

Seismic velocity
Seismic velocity was a parameter included in several processing steps (e.g., post-STM, depth conversion, and thermal conductivity estimation).To assess the uncertainty of heat flow derived from seismic velocity, we validated the interval velocity obtained by velocity conversion.Since uncertainty in the P-wave velocity within the sequence above the BSR may result in apparent BSR depths that too shallow or too deep (Ohde et al. 2018), we evaluated P-wave velocity within the shallow sequence.We first presumed that the P-wave velocity range in methane-hydrate-bearing sediment was 1800-2300 m/s on the basis of our results of seismic velocity analysis, and then we calculated the BSR depth and the heat flow using velocities of 1800 m/s and 2300 m/s within the corresponding depth in the time domain (two-way-travel time of 0.5 s below the seafloor).Given that we adopted a constant value of thermal conductivity, this velocity variation corresponds directly to a 20% uncertainty in calculated heat flows.

Assumed thermal conductivity
Since we based our constant value of thermal conductivity (1.2 W/mK) on the empirical relationship between P-wave velocity and thermal conductivity (Duffaut et al. 2018), this average value of thermal conductivity differs from the high-resolution record derived from borehole data (Ganguly et al. 2000).To evaluate the resulting uncertainty, we referred to measured thermal conductivity data from Ocean Drilling Program Sites 808 (Fisher 1993) and International Ocean Drilling Program Sites C0002 (Sugihara et al. 2014).Since the measured thermal conductivity range in such sediment was 1.0-1.5 W/mK, the uncertainty in thermal conductivity may be as great as 25% and would lead to a similar uncertainty of heat flow.

Depth error derived from migration processing approaches
Although we performed post-STM in the interest of saving computation costs, seismic profiles obtained from pre-stack time migration (pre-STM) are more accurate than those obtained from post-STM.We evaluated the difference in heat flow derived from the two methods in the shallow sequence above the BSR.
As shown in Fig. 7a and b, we created seismic profiles using post-STM and pre-STM, and then we converted them to the depth domain using the same P-wave velocity model, which is produced by the high-density velocity analysis noted above, to evaluate the error due to the migration processing approach.Finally, we calculated the heat flow in each profile, which is compared in Fig. 7c.Even though the resultant depth profiles agree with the seafloor depth, the BSR depth differed slightly in several parts of the profile (Fig. 7d), corresponding to a difference in heat flow as great as 13%.We concluded that the uncertainty due to the migration processing approach reaches ± 13%.
Although the total uncertainty integrated individual uncertainty is higher than those of similar studies (Ganguly et al. 2000;Townend 1999;Ohde et al. 2018), the   5a, respectively.The yellow translucent area with depth contours indicates the distribution of plutonic rocks, inferred from on-land exposures and apparent beneath an unconformity between the cover sediment and the igneous basement on the seismic profiles (Kimura et al. 2022) effects of some errors cancel and relative errors comparing profile to profile is perhaps less than 10% (Marcaillou et al. 2006).

Factors of high heat flow anomalies
Although the map of BSR temperature (Fig. 6a) shows moderate values and few if any apparent anomalies, the map of heat flow (Fig. 6b) shows notable anomalies on each seismic profile.The reasons for these have been discussed in previous studies (e.g., Martin et al. 2004;Yamano et al. 1992;Ohde et al. 2018).
The inferred high heat flow anomalies shown in Figs.5a  and b appear to be associated with seafloor erosion and fluid flow heat advection along active faults (e.g., Martin et al. 2004;Kunath et al. 2021).In Fig. 5a, the BSR in the seaward portion of KIN_054 is not parallel to the seafloor and is likely modified by surface erosion along the seaward flank of the prominent anticline (Fig. 6d).Since the BSR depth is shallow, heat flows in the vicinity of surface erosion have high values.In contrast, in Fig. 5b, the sediment suppresses heat flow variations.On the basis of our results of seismic profile and heat flow estimates, we concluded that the spatial variation of heat flow implies anomalous characteristics in the subsurface.Here, we propose two influencing factors that are consistent with the presence of active faults: erosion and fluid migration.Martin et al. (2004) pointed out that heat flow anomalies derived from BSR can be caused by erosion and sedimentation.Our results indicate extremely high heat flow values, approximately 100 mW/m 2 , in the accretionary prism (Fig. 4d).Such anomalies can be caused by erosion and could rival or exceed those caused by fluid flow (Martin et al. 2004) (Fig. 5).Despite our inability to confidently apportion between erosion and fluid flow as causes for the observed heat flow anomalies, at least some erosion-induced heat flow anomaly is expected to exist here because this tectonic setting is known to contain active faulting, which has shown to enhance erosion along the seaward flanks of growing anticline seafloor ridges (Alves et al. 2014).

Erosion
A general assumption that the BSR marks a stable thermal condition and an equivalence to the BGHS (Grevemeyer and Villinger 2001; Yamada et al. 2014) can be invalid where erosion or sedimentation has occurred (Martin et al. 2004).Tectonic factors which generate topographic uplift can thus produce non-equilibrium temperature profiles that yield uncertainties in heat flow (Mandal et al. 2014).After a disturbance, the BGHS migrates upward or downward in response to any changes in temperature and pressure conditions.Delisle et al. (1998) reported BSR downward migration due to rapid erosion (slope failure) and upward migration due to rapid sedimentation, resulting in errors in heat flow values of roughly 5%.Kinoshita et al. (2011) have reported that the time for the BSR to re-equilibrate after such a disturbance is ~ 10 kyr in the Nankai accretionary prism.Therefore, the BSR temperature profiles in regions affected by tectonic disturbances might be subject to non-equilibrium conditions that affect the heat flow value.
In addition, several factors (e.g., tectonic uplift and gas sources) can affect the depth of BGHS and thus constrain the BSR positions, including double BSRs (Matsumoto et al. 2004;Zhang et al. 2022).In our study area, we interpret instances of recent BSR downward migrations near sites of erosion (Fig. 8) as well as eroded areas without a corresponding migration and thereby can interpret erosion as the cause of the BSR migration.As mentioned above, BSR migration can take place for several thousands of years to keep thermal equilibrium conditions again.Thus, if the eroded areas without a corresponding migration reflect non-equilibrium conditions, then the extremely high heat flow anomalies observed may be consistent with recent or ongoing fault activity.

Fluid migration
We detected local high heat flow anomalies in the ITZ at locations where erosion cannot be unambiguously identified, and for this reason, we considered that BSR depth may be perturbed by fluid flow causing heat advection (Fig. 5b, d).Fluid flow in the accretionary prism could partially originate in the décollement and rise toward the seafloor into the ITZ, where a BSR forms (Yamano et al. 1992;Baba and Yamada 2004).Our velocity models show low-velocity zones accompanying frontal or out-of-sequence thrusts (Fig. 3).This observation is consistent with fluid flow along major thrust faults (Fig. 9).Moreover, in the region of pronounced heat flow anomalies, thrust faults appear to be exposed at the seafloor, whereas in areas with low heat flows, the faults are blanketed by slope sediments.Since sediment cover is likely to impede or prevent fluid flow (Kobayashi 2002), slope sediments in the ITZ may act to impede fluid flow through the thrust faults.Since our results of heat flow reflect the BSR depth, the area without high heat flow anomalies may indicate an absence of fluid flow along a fault.Although a conclusive determination of the activity status of sediment-covered faults cannot be reliably made based solely on heat flow analysis, the conspicuous absence of cover sediment deformation resulting from fault displacement suggests that the faults in this area could be inactive (Fig. 9) or sedimentation is very rapid.
We classified high heat flow anomalies on the basis of heat flow values, the nearby topography, and the seismic profiles.Where the seafloor appeared to be eroded above faults near a high heat flow anomaly, we assigned it an erosional origin.This type of anomaly also tended to have higher heat flows than anomalies suspected to have been caused by fluid migration.High heat flow anomalies attributed to fluid flow tended to be in the accretionary prism and to be surrounded by low heat flow anomalies, where the presence of slope sediments may inhibit fluid flow.In several areas, we cannot estimate the heat flow, nor can we infer or rule out the presence of active faults due to the absence of the BSR, which might be associated with tectonic activities such as surface erosions.Nevertheless, the estimated locations of active faults are consistent with the results of previous studies (e.g., Tsuji et al. 2014).

Landward heat flow trends
Our heat flow distribution shows a general trend of decreasing from ~ 50 mW/m 2 to < 40 mW/m 2 landward (Fig. 6b), which reflects the processes of thickening accretionary wedge sediments and subduction of the Philippine Sea plate.This trend is consistent with previous studies (e.g., Ashi et al. 2002b;Hamamoto et al. 2011;Harris et al. 2013).One noteworthy feature is that the heat flow on Line KIN_030, the westernmost profile, differs from the other survey lines in having some exceptionally low heat flow values.A subducted seamount has been detected beneath the accretionary prism to the west of this profile (e.g., Park et al. 1999;Kodaira et al. 2000;Gulick et al. 2004;Nakamura et al. 2022), which may affect the heat flow profile in two contrasting ways.Spinelli and Harris (2011) have pointed out that the subduction of a seamount would impede fluid flow by closing original pathways and consequently lessen the advective transfer of heat.Likewise, the closing of pathways by the seamount subduction may suppress the migration of heated fluid from greater depth to the shallower part in this area and result in a redistribution of heat that would make the landward part lower heat flow with main contribution of conduction in the accreted sediment.In contrast, other studies have argued that seamount subduction may cause fracturing of the upper plate (e.g., Dominguez et al. 1998;Wang and Bilek 2011), in which case these fractures might act as fluid paths and enhance heat flow.Since our seismic profile did not image pervasive fracturing, we interpret the observed low heat flow on profile KIN_030 to be the result of reduced advective heat transfer.
As mentioned above, heat flow generally decreases landward; nevertheless, our results show a zone of relatively anomalous heat flow (> 40 mW/m 2 ) in the eastern Near the eroded seafloor, there may be two BSRs: one that represents the pre-erosional pressure and thermal regime (black dashed line at the rightmost part), and another that has since diverged away from the seafloor in response to the post-erosional pressure and thermal regime (black line).Such surface disturbances can also contribute to high heat flow anomalies part of the study area (Fig. 6b, d).The extent of this trend in our study area seems to correspond with the presence of plutonic rocks beneath the Kii Peninsula, which may extend to the offshore (Kodaira et al. 2006;Kimura et al. 2014Kimura et al. , 2022;;Arnulf et al. 2022).Since the thermal conductivity of plutonic rocks is generally higher than that of sedimentary rocks (Clauser and Huenges 1995) and the regional heat flow is derived from the subducting Philippine Sea plate beneath the onshore and offshore portions of the Kii Peninsula (e.g., Ji et al. 2016), the thermal structure in this region can be constrained by the existence of the plutonic rock.Therefore, the higher thermal conductivity of the plutonic rock may contribute to the high heat flow trend observed in the landward side in the eastern part.The Kii Peninsula contains hot springs with high temperatures, high geothermal gradients, and heat flow anomalies even though this area is not a volcanic region (Furukawa et al. 1998;Yano 1999;Tanaka et al. 2004).These geothermal anomalies may be associated with extensional faults and vertical dykes with N-S strikes (Umeda et al. 2006).Even though there may be still uncertainty of the heat flow profile in our study area, we concluded that the presence of the plutonic rock is contributing to high heat flow trend in our study area.

Conclusions
We conducted seismic processing with automated velocity analysis of parallel seismic profiles acquired during the 2018 R/V Kaimei voyage (KM18-10) offshore the Kii Peninsula.Using high-resolution P-wave velocity information, we mapped the distribution of the BSR in the study area.Velocities in the methane hydrate-bearing sediments above the BSR ranged from 1800 to 2300 m/s and allowed us to outline low-velocity zones along the major thrust faults that suggest the presence of advective fluid flow.Heat flow, calculated from the BSR depth, ranged from 20 to 100 mW/m 2 , although this range includes uncertainties due to seismic velocity, migration processing approach, and thermal conductivity.Our heat flow distribution map shows high and low anomalies, which can be attributed to seafloor erosion or fluid flow.
High heat flow anomalies due to erosion, which can be confirmed by seismic profiles or bathymetric mapping, can exceed those due to fluid flow.These anomalies can serve as evidence of recent tectonic activity that confirms the existence of active faults.Low heat flow anomalies, detected beneath areas of thick cover sediment on our seismic profiles, are interpreted as the consequence of rapid sedimentation.Areas of high heat flow neighboring these low anomalies are usually not covered by slope sediments and are interpreted as the result of heat advection by fluid flow through thrust faults that reach the seafloor.Since active faults can offer highly permeable fluid paths, active faults may be associated with high heat flow anomalies caused by fluid flow.We also found evidence that parts of the BSR in the study area may represent conditions out of thermal equilibrium.
We noted a trend toward higher heat flow on the landward part of our study area close to the Kii Peninsula.We interpret this higher heat flow anomaly results from the thermal structure of plutonic rocks beneath the Kii Peninsula, which may extend offshore.

Fig. 1
Fig. 1 Bathymetric map of the landward side of the Nankai Trough off the Kii Peninsula.The black triangles along the plate boundaries on the upper panel denote the subduction trend.R/V Kaimei KM18-10 seismic survey lines are shown in red.The black dashed region is the location of subducted seamount inferred by Kodaira et al. 2000.White dashed lines are the segment of megathrust rupture along the Nankai Trough(Ando 1975;Cummins et al. 2002).The white star indicates the location of Ocean Drilling Program Sites 808(Mikada et al. 2002)

Fig. 3
Fig. 3 P-wave velocity structure for all six survey lines in this study.(a) Line KIN_030.(b) Line KIN_038.(c) Line KIN_046.(d) Line KIN_54.(e) Line KIN_062.(f) Line KIN_070.The area with weak reflections is masked in gray.The expanded views on (c) and (d) show examples of low-velocity zones.The red and black arrows indicate anomalous low-velocity zones corresponding to thrust faults and relatively low-velocity area below BSRs, which may reflect the presence of the free gas, respectively.The white arrowheads indicate BSR positions

Fig. 4
Fig. 4 Seismic depth profile (above) and semilogarithmic heat flow plot (below) for the six survey lines in this study.(a) Line KIN_030.(b) Line KIN_038.(c) Line KIN_046.(d) Line KIN_54.(e) Line KIN_062.(f) Line KIN_070.Blue dashed rectangles on (f) correspond to the regions shown in Figs. 5 and 8.The Imbricate Thrust Zone (ITZ) was formed in the middle of the accretionary prisms.The plutonic rock (PR) inferred by Kimura et al. 2022 may be underlain in (e) and (f).Yellow arrowheads denote the BSRs.Red arrows indicate high heat flow anomalies associated with eroded seafloor

Fig. 5
Fig. 5 (a), (b) Detail of Line KIN_054 seismic profile and (c), (d) corresponding heat flow plot; location in Fig. 4. Yellow arrowheads denote the BSR.Red arrowheads in (b) indicate the relic and migrated BSRs caused by erosion.Red lines and orange translucent zone in (a) and (b) indicate the interpreted faults and cover sediment, respectively.Red and blue arrows in (c) and (d) denote local increasing and decreasing trends of heat flow, respectively.Gaps in BSR-derived heat flow data exist in (d) because heat flow cannot be estimated where the BSR is absent or weak.Note that the high heat flow anomaly around 50 km in (b) corresponds to a feature in the seafloor topography

Fig. 6
Fig. 6 Distribution of (a) BSR temperature and (b) heat flow off the Kii Peninsula.Thick and thin lines show our results and the results of Ohde et al. (2018), respectively.Note that the range of the color scale in (b) is limited to 30-60 mW/m 2 so as to emphasize heat flow anomalies caused by fluid migration.Higher values are omitted.(c) Location of active faults related to erosion or fluid migration, which are inferred from our interpretation of topography heat flow, and seismic profiles.The topographic feature is shown as slope angles derived from high-resolution bathymetry.(d) Heat flow in this study area with geological features superimposed.The black dashed circle and yellow dashed region indicate the zone with the anomalous high heat flow and the extremely high heat flow area (EHFA), shown in Fig.5a, respectively.The yellow translucent area with depth contours indicates the distribution of plutonic rocks, inferred from on-land exposures and apparent beneath an unconformity between the cover sediment and the igneous basement on the seismic profiles(Kimura et al. 2022)

Fig. 7
Fig. 7 Extracted seismic profiles obtained by (a) pre-STM and (b) post-STM.(c) The difference in heat flow between (a) and (b).(d) The difference between BSR depth and seafloor depth in profiles obtained by pre-STM (black) and post-STM (red).Note the different vertical scales of (c) and (d).Yellow arrowheads in (a) and (b) represent the location of BSRs.The main factors for the uncertainty can be caused by the difference of the BSR depths

Fig. 8
Fig. 8 Detail of the seismic profile of Line KIN_054.Yellow arrowheads indicate the BSR.The orange dashed and solid regions indicate the BSR before and after depth change constrained by thermal conditions, respectively.The edge of the BSR has moved downward at its seaward end in response to plausible surface erosion, though it also could be caused by other factors

Fig. 9
Fig. 9 Schematic image of geological features and heat flow distribution.The orange arrows represent the fault displacement.Where slope sediment covers thrust faults, low heat flow is found.At such sites, the sediment cover might prevent fluid flow (small blue arrows).Areas without cover sediment and with eroded surfaces have high heat flow.Around these areas, fluid flow may reach the seafloor through the active faults (large blue arrows).Near the eroded seafloor, there may be two BSRs: one that represents the pre-erosional pressure and thermal regime (black dashed line at the rightmost part), and another that has since diverged away from the seafloor in response to the post-erosional pressure and thermal regime (black line).Such surface disturbances can also contribute to high heat flow anomalies