Surface heat flow measurements in the eastern margin of the Japan Sea using a 15 m long geothermal probe to overcome large bottom-water temperature fluctuations

Accurate surface heat flow data are required for a wide range of geological and geophysical applications. However, sediment temperature measurements beneath the seafloor often involve large uncertainties owing to the influence of bottom-water temperature (BWT) fluctuations. Previous studies reported apparently negative geothermal gradients in the Joetsu Basin of the Japan Sea and suggested that BWT fluctuations disturbed sediment temperatures. To address this problem, we monitored BWTs in the Joetsu Basin over a 2 year period to determine the depth at which the influence of BWT fluctuations on sediment temperature becomes negligible. Combined with sediment thermal diffusivity data, we determined that the BWT fluctuations can disturb sediment temperatures to a depth of 2 m. We obtained heat flow values of 81–88 mW m− 2 by measuring sediment temperatures at depths > 2 m using a 15 m long geothermal probe. The measured heat flow values are inversely correlated with topography owing to the effect of topographic change on the geothermal structure near the seafloor. A two-dimensional geothermal structure model was constructed to account for the topography, yielding an estimated regional background heat flow of 85 ± 6 mW m− 2. This study provides two important guidelines for obtaining accurate surface heat flow data in marine areas with large-amplitude BWT fluctuations: (1) quantitative information regarding BWT fluctuations and sediment thermal diffusivity is required to evaluate the depth range to which BWT fluctuations affect sediment temperature; and (2) information regarding the lithology and consolidation state of seafloor sediments is required for effective penetration using a long probe.


Introduction
The flow of heat from the Earth's deep interior to its surface can be used to constrain subsurface geothermal structures (e.g., Hamamoto et al. 2011;Harris et al. 2017;Suenaga et al. 2018) and fluid flow regimes (e.g., Fisher et al. 2005;Hutnak et al. 2006;Harris et al. 2020). Heat flow is essentially the amount of thermal energy that passes through a unit area over a unit time, and is calculated as the product of the vertical geothermal gradient and thermal conductivity of the material over the depth interval where the geothermal gradient is obtained. Geothermal gradients in marine areas are generally calculated from sediment temperatures obtained using a geothermal probe (typically 3-5 m in length) with several sensors to measure temperatures at different depths. The thermal conductivity of sediment is measured in the laboratory on sediment core samples (Von Herzen and Maxwell 1959;Sass et al. 1984;Goto et al. 2017) or in situ (Hyndman et al. 1979;Villinger and Davis 1987;Hartmann and Villinger 2002).
A critical aspect to obtaining reliable heat flow data in marine areas is that the sediment temperature measurements must be made under stable bottom-water temperature (BWT) conditions, such as in deep-sea areas. Seafloor sediment temperatures temporally change in regions with large BWT fluctuation amplitudes (e.g., shallow-sea areas), which undermine the accuracy of heat flow measurements (Davis et al. 2003). For example, numerous marine heat flow data have been obtained in deep-sea areas around Japan, whereas heat flow data in shallow-sea areas near the coast remain limited (Tanaka et al. 2004). To evaluate the effect of BWT fluctuations on the geothermal gradient, Yamano et al. (2014) monitored BWTs for more than 1 year at three stations in the Nankai accretionary prism southeast of Kii Peninsula, southwest Japan. They indicated that BWT fluctuations with a maximum temperature change of 0.39 K can produce geothermal gradient errors in the order of 20-30% using a 3 m-long geothermal probe in areas with water depths of approximately 2500-2800 m. The effect of BWT fluctuations on the geothermal gradient measurements obtained using the same probe was reduced to 10% at ~ 3300 m water depth, where a maximum temperature change of 0.1 K was observed.
As BWT fluctuations propagate through sediment via heat conduction, their amplitude decreases exponentially with increasing depth below the seafloor (e.g., Carslaw and Jaeger 1959). The penetration depth of BWT fluctuations in sediment depends on the predominant fluctuation period and thermal diffusivity of the sediment. Longerperiod BWT fluctuations penetrate deeper into sediment, and BWT fluctuations propagate deeper into sediment with higher thermal diffusivity . To obtain reliable heat flow measurements in marine areas where the sediment temperature is affected by BWT fluctuations, it is therefore necessary to either (i) measure sediment temperatures at different depths where the effect of BWT fluctuations has sufficiently attenuated or (ii) remove the effect of BWT fluctuations from the sediment temperature data obtained at different depths. In the former case, the sediment temperatures undisturbed by BWT fluctuations can be measured in a borehole (e.g., Harris et al. 2011). However, this method is generally difficult to apply owing to high drilling costs. In the latter case, sediment temperatures undisturbed by BWT fluctuations can be obtained by monitoring both bottom-water and sediment temperatures over a long period and by removing the BWT fluctuation effect from the sediment temperatures. Although this method is useful for obtaining reliable heat flow data, a special instrument is required to simultaneously measure bottom-water and sediment temperatures over a long period of time, such as that developed by Kinoshita et al. (1996) and Hamamoto et al. (2005), which limits the number of heat flow measurements. Heat flow can also be indirectly estimated from seismic profiles that show bottomsimulating reflectors (BSRs) associated with gas hydrates (Yamano et al. 1982;Hyndman et al. 1992). Although this approach is useful for estimating heat flow over a wide area where BSRs are observed on a seismic profile, the estimated heat flow has a relatively large uncertainty (Grevemeyer and Villinger 2001).
In this study, we propose a method to obtain heat flow at the seafloor (hereafter, surface heat flow) that is undisturbed by BWT fluctuations, by measuring sediment temperatures using a long geothermal probe that is lowered from a research vessel to deep sub-bottom depths where the effect of BWT fluctuations has sufficiently decayed. The subbottom depth at which sediment temperatures are affected by BWT fluctuations is evaluated using a time-series of BWT data obtained in the area. Surface heat flow is calculated from the sediment temperatures measured at depth where the effect of BWT fluctuations is sufficiently negligible. The application of multi-penetration (i.e., pogo-style) measurements permits reliable surface heat flow data to be efficiently obtained. This method was applied in the Joetsu Basin ( Fig. 1) in the eastern margin of the Japan Sea where negative geothermal gradients were reported (Machiyama et al. 2009), likely owing to BWT fluctuations, to investigate the heat flow supplied from substantial depth (i.e., regional background heat flow) and to infer the geothermal structure of the basin. Sediment temperatures were measured using a 15 m long geothermal probe. Bottom-water temperatures were measured over a 2 year period to evaluate the depth to which BWT fluctuations disturb the sediment temperature, and surface heat flow values undisturbed by BWT fluctuations were successfully obtained.

3
We first present the results of the long-term BWT monitoring study and evaluation of the effect of BWT fluctuations on the temperature of sediment near the seafloor in the Joetsu Basin. We then show the surface heat flow results determined from undisturbed sediment temperature data obtained using a 15 m long geothermal probe. The obtained surface heat flow values seem to inversely correlate with topography owing to topographic effects on the geothermal structure near the seafloor (e.g., Ganguly et al. 2000;He et al. 2007). We estimate a range of regional background heat flow values supplied from deep in the basin by creating a two-dimensional (2D) geothermal structure model that considers topographic change using the measured seafloor heat flow values as thermal constraints at the seafloor. We present two important guidelines for obtaining accurate surface heat flow data in marine areas with large-amplitude BWT fluctuations using a long geothermal probe.
In the present study, we demonstrated that reliable surface heat flow data undisturbed by BWT fluctuations are obtained by penetrating a long geothermal probe into sediment in the Joetsu Basin where large-amplitude BWT fluctuations were observed. This method can be used for efficient measurement of heat flow in marine areas where BWT fluctuations prevent reliable measurement of heat flow, such as shallow sea areas near the coast.

Geologic setting
The Japan Sea is a semi-isolated, back-arc basin surrounded by the Japan island arc and Eurasian continent (Fig. 1). The Japan Sea opened between 32 and 10 Ma (Tamaki et al. 1992;Jolivet et al. 1994), at which time subsidence occurred in the sea and its margins (Ingle 1992). The basin is filled with thick sediments including good petroleum source rocks (Suzuki 1989;Muramoto et al. 2007). The eastern margin of the Japan Sea has been subjected to E-W compressional stress since the Pleistocene (Tamaki et al. 1992;Jolivet et al. 1994), and a series of NE-to SW-trending anticline-syncline systems have formed along its eastern margin (Okamura et al. 1995).
The Joetsu Basin ( Fig. 1) is located to the west of a major oil and gas field along the coast of the Japan Sea in northeastern Japan and is expected to be potentially rich in fuel resources (Osawa et al. 2005;Monzawa et al. 2006Monzawa et al. , 2011Okui et al. 2008;Nguyen et al. 2016). Six major formations of mainly mudstone and sandstone have been identified above the Green Tuff, supposedly from the regional volcanic basement in the basin: the Nanatani Formation (Early to Middle Miocene); Lower Teradomari Formation (Middle Miocene); Upper Teradomari Formation (Late Miocene); Shiiya Formation (Early Pliocene); Nishiyama Formation (Late Pliocene to Pleistocene); and Haizume Formation (Pleistocene) (Muramoto et al. 2007;Okui et al. 2008). The  Kisimoto (2000). Rectangle shows the study area. Detailed bathymetry of the area is shown in Fig. 2. The thick solid line crossing the study area indicates part of the METI 2D seismic survey line (Japan National Oil Corporation 2003). Index map is inserted 1 3 basin is characterized by two major topographic features formed by asymmetric anticline structures (Fig. 2). One is a spur that sticks out from the continental slope, called the Umitaka Spur; the other is an isolated NE-SW-trending knoll, named the Joetsu Knoll, located 8-10 km northwest of the Umitaka Spur. Several mounds (diameter ~ 450 m; relative elevation ~ 30 m) and pockmarks (diameter ~ 600 m; relative depth ~ 50 m) occur on the summit areas of the Umitaka Spur and Joetsu Knoll (Matsumoto et al. 2009;Saeki et al. 2009;Freire et al. 2011). Seismic profiles across the Umitaka Spur and Joetsu Knoll reveal weak or blanking reflections below the mounds and pockmarks, indicating upward fluid migration controlled by the anticlinal fault system (Muramoto et al. 2007;Matsumoto et al. 2009;Saeki et al. 2009;Freire et al. 2011;Nakajima et al. 2014). BSRs have been identified on seismic profiles obtained from the Umitaka Spur and Joetsu Knoll (Matsumoto et al. 2009;Saeki et al. 2009;Freire et al. 2011;Nakajima et al. 2014). Gas hydrates in mud-dominant sediment have been recovered in and around the mounds and pockmarks (Matsumoto et al. 2009;Hachikubo et al. 2015). Active cold seep sites with biological communities, outcrops of massive gas hydrates, and active methane gas seeps have been observed on several mounds on the Umitaka Spur and Joetsu Knoll (Matsumoto et al. 2009). Machiyama et al. (2009) intensively measured surface heat flow in the Joetsu Basin to investigate the pore fluid flow in sediment associated with cold seep activity in the basin. Prior to heat flow measurements by Machiyama et al. (2009), heat flow data were unavailable in the Joetsu Basin. Machiyama et al. (2009) used two types of geothermal probes to measure the geothermal gradient: (i) a piston corer with several small temperature loggers mounted on the barrel of several meters in length; and (ii) a stand-alone heat flow meter (SAHF) with five thermistors inside a probe of 60 cm in length (Kinoshita et al. 2006). The SAHF was developed to measure sediment temperatures using a submersible. Machiyama et al. (2009) used these instruments and obtained surface heat flow values ranging from 65 to 106 mW m -2 in areas without active cold seep sites in the Joetsu Basin, and high surface heat flow values of > 150 mW m -2 were obtained at the active cold seep sites on the Umitaka Spur and Joetsu Knoll. Machiyama et al. (2009) suggested that advective heat transport owing to upward fluid flow contributed to produce the high surface heat flow values. On the contrary, negative temperature profiles were observed at and around the active cold seep sites on the Umitaka Spur and Joetsu Knoll using the SAHF. Machiyama et al. (2009) suggested that the negative temperature profiles were caused by large-amplitude BWT fluctuations.

Sediment temperature disturbances due to BWT fluctuations
To monitor the BWT fluctuations, a water temperature recorder NWT-DN (NiGK Corporation, Kawagoe, Japan) was deployed at the summit area of the Joetsu Knoll (Station WT in Fig. 2) on June 17, 2010 and recovered on May 26, 2012 using a remotely operated vehicle. The accuracy and resolution of the logger are 0.05 and 0.001 K, respectively. The sampling interval was set to 20 min. 2 years of BWT data were thus obtained (Fig. 3). The recorded BWT values show long-period fluctuations with sudden large-amplitude changes (~ 0.11 K), likely reflecting the replacement of bottom-water mass. Such long-period and large-amplitude fluctuations are characteristic for conditions that tend to substantially disturb the sediment thermal field below the seafloor.
We calculated the sediment temperature variations expected from the observed BWT fluctuations to evaluate the effect of BWT fluctuations on the sediment thermal field in the Joetsu Basin. We assumed semi-infinite sediment with uniform thermal properties. As a boundary temperature condition, we applied a series of step functions that could approximate non-periodic, sudden, large-amplitude BWT changes. Sediment temperature is expressed by the following analytical solution (Carslaw and Jaeger 1959;Goto and Yamano 2010): (1) where T(z, t M ) is the sediment temperature at depth z below the seafloor at time t M , T 0 is the average BWT, G is the geothermal gradient, ΔT i is the temperature change relative to T 0 in the time between t i−1 and t i prior to t M , κ s is the thermal diffusivity of the sediment, and erfc(x) is the complimentary error function defined as: To evaluate the effect of BWT fluctuations on the sediment temperature, we assumed G = 100 mK m -1 and κ s = 2.39 × 10 -7 m 2 s -1 , which is the average thermal diffusivity of sediment at a depth of 0-5 m in the Joetsu Basin (Goto et al. 2017). Figure 4 shows synthetic sediment temperatures at various depths calculated from the assumed G and κ s values alongside the observed BWT fluctuations. The sediment temperatures reflect values after the BWT fluctuation effect arrived at those depths. Figure 5 shows synthetic sediment temperature profiles at various times. These figures demonstrate significant thermal disturbances of up to 0.02 K owing to the BWT fluctuations at depths of 1.5 m below the seafloor. However, the thermal disturbance is less than 0.02 K for synthetic temperatures deeper that 2 m below the seafloor. The modeling shows a potential error range of − 63% to + 75% for the geothermal gradient calculated from sediment temperatures down to 60 cm below the seafloor, − 19% to + 23% for the gradient calculated from temperatures down to 2 m below the seafloor, ± 7% for the gradient calculated from temperatures down to 5 m below the seafloor, and ± 2% for the gradient calculated from temperatures in the depth range of 2-5 m below the seafloor. The results indicate that sediment temperatures must be measured at several points at a depth greater than 2 m to obtain a geothermal gradient without thermal disturbances owing to BWT fluctuations in the Joetsu Basin.
(2) Fig. 4 Synthetic sediment temperature variations expected from BWT fluctuations recorded at station WT on the Joetsu Knoll (Fig. 3). An undisturbed geothermal gradient of 100 mK m -1 and sediment thermal diffusivity of 2.39 × 10 − 7 m 2 s -1 were assumed for the calculations 1 3

Surface heat flow measurements
To investigate the distribution and formation mechanism of gas hydrates in mud-dominant sediments at relatively shallow depths below the seafloor, the MD179 cruise was carried out in June 2010 using the French research vessel Marion Dufresne in the Joetsu Basin. Surface heat flow measurements were conducted at seven stations aligned along the METI 2D seismic survey line (Japan National Oil Corporation 2003), which crosses the southwestern part of the Joetsu Knoll and northern to northeastern part of the Umitaka Spur (Fig. 2). Several cold seep activities were observed on the seafloor of the Joetsu Knoll and Umitaka Spur (Matsumoto et al. 2009). Fluid migrations in and around such cold seep sites should affect local geothermal structure. To avoid possible thermal disturbance by migrating fluid, we carried out surface heat flow measurements in areas without cold seep activity.
Sediment temperatures for calculating surface heat flow were measured using a giant sediment core sampler of the R/V Marion Dufresne, as shown in Fig. 6. A 15 m long geothermal probe was prepared by welding a conical steel to one side of a core barrel; thus no sediment was recovered by the penetration. Seven miniaturized temperature data loggers (ANTARES Datensysteme GmbH, Stuhr, Germany; see Pfender and Villinger (2002) for a detailed description) were attached to the probe at an interval of 2 m. The geothermal probe was connected to a weight head of 4.5 or 2.5 t for penetration into the sediment. The tilt was measured during the probe penetration to correct the sensor depth using a tilt meter (Model FL-300-0002, Kaiyo Denshi Co., Ltd., Tsurugashima, Japan). During the MD179 cruise, two probe penetrations were conducted in a pogo-style fashion at each station. The BWT values were measured using the temperature data loggers of the probe prior to penetration into the sediment. Frictional heat due to the probe penetration into the sediment was corrected using the linear regression of temperature versus the inverse time elapsed since the probe penetration (Hyndman et al. 1979). We obtained sediment temperature profiles with a geothermal gradient in the range of 88-97 mK m -1 ( Fig. 7; Table 1). Unfortunately, temperature data were not obtained from the two upper data loggers at stations HF1 and HF2.
We estimated sub-bottom depths of the individual temperature loggers assuming that (i) the BWT values were stable and (ii) the average BWT value measured prior to probe penetration was equal to the temperature on the seafloor. The estimation procedures were as follows. We first linearly extrapolated the sediment temperature profile to the seafloor. We then found the topmost temperature logger depth where the extrapolated temperature equaled the average measured BWT. The sub-bottom depths of the other temperature loggers were then calculated. The estimated sensor depths were greater than 2.3 m below the seafloor (Fig. 7). The maximum depth of the sediment temperature measurements was estimated to be 16.3 m below the seafloor, indicating that the penetration of the probe into the sediment reached the lower part of its weight head because of its heavy weight. However, the BWT values obtained in the Joetsu Basin fluctuated (Fig. 3), indicating inconsistency with assumption (i) that the BWTs were stable. We evaluated the uncertainties of the sensor depths using the standard deviation of the BWT fluctuation (0.018 K) as ~ 0.20 m by applying the same estimation procedure explained above. As shown in Fig. 4, sediment temperatures at depths less than 2 m below the seafloor are expected to be disturbed by BWT fluctuations in the Joetsu Basin. By considering the uncertainty of the estimated sensor depth, all of the sediment temperatures in Fig. 7 can be used to calculate the surface heat flow undisturbed by BWT fluctuations.
The Joetsu Basin may have longer-period BWTs that cannot be detected from the 2 years of BWT monitoring. If there are such longer-period BWT fluctuations with large amplitudes in the basin, sediment temperatures at depths greater than 2 m below the seafloor will be affected by the fluctuations, and sediment temperatures measured with a long geothermal probe will show a non-linear profile. Sediment temperatures measured with a 15 m long geothermal probe in the present study show almost linear profiles (Fig. 7), suggesting that there was no or little effect of longer-period BWT fluctuations on the sediment temperatures.
Thermal properties, such as thermal conductivity, heat capacity, thermal diffusivity, and index properties, such as porosity, density, of the sediment were measured from core samples recovered at eight stations (Fig. 2) using the R/V Marion Dufresne corer during the MD179 cruise (see Goto et al. (2017) for more details). The obtained thermal conductivity data showed large variability. The thermal conductivity of the sediment near the seafloor was found to be approximately 0.8 W m -1 K -1 , and to increase with depth, reaching ~ 1.0 W m -1 K -1 at a depth of 30 m below the seafloor. Porosity decreases exponentially with depth over the same depth interval, whereas bulk density increases. The observed increase in thermal conductivity and bulk density can be explained by compaction owing to sediment accumulation (Goto et al. 2017).
The sediment thermal conductivity data by Goto et al. (2017) were used to calculate surface heat flow. The Bullard plot method (Bullard 1939) was applied to account for the change of thermal conductivity with depth in the surface heat flow calculation (see Appendix for calculation details). The calculated surface heat flow values range from 81 to 88 mW m -2 (Table 1). Two heat flow values obtained at each station  The distribution of surface heat flow along the METI 2D seismic profile is represented in Fig. 8. Except for station HF1, the surface heat flow values seem to inversely correlate with topography; surface heat flow is relatively high near the base of the slope of the Joetsu Knoll and Umitaka Spur and relatively low on the Joetsu Knoll and Umitaka Spur. Because of topographic effects on the subsurface thermal structure, the conductive heat flow focuses (defocuses) over concave (convex) topographic regions, resulting in apparently higher (lower) surface heat flow values (Lachenbruch 1968;Blackwell et al. 1980;He et al. 2014). The surface heat flow variation (except for station HF1, which is discussed below) can therefore be attributed to topographic changes. Station HF1 is located at the base of a steep scarp (Fig. 2) and its surface heat flow value is expected to be relatively high owing to the topographic effect on heat flow on the seafloor; nevertheless, low surface heat flow values were obtained ( Fig. 8; Table 1). The seismic profile in Fig. 8 suggests a sediment supply from the scarp to the seafloor around station HF1. Rapid sediment accumulation at the seafloor is known to reduce surface heat flow (Hutchison 1985;Wang and Davis 1992). Such sedimentation effects are therefore likely to be the cause of the low surface heat flow values at station HF1.  (i) at least three sediment temperatures were measured at depths greater than 2 m below the seafloor; and (ii) the temperature profile, including temperatures measured at depths shallower than 2 m below the seafloor, is linear. According to these criteria, all of the temperatures measured using a SAHF were likely disturbed by BWT fluctuations because of the short probe length (60 cm). Surface heat flow data obtained using a SAHF were thus not used in this study. For the geothermal gradient data measured by Machiyama et al. (2009) using a piston corer, eight measurements collected on and around the Umitaka Spur satisfied the criteria (Fig. 2). In these measurements, the piston corer barrel penetrated the sediment to a depth of ~ 5 m below the seafloor. These surface heat flow values are equal to or are greater than those measured in the present study. A high surface heat flow value (101 mW m -2 ) was obtained at the base of the steep western slope of the Umitaka Spur (Fig. 2), which likely reflects the topographic effect. Another high surface heat flow value in the basin floor west of the Umitaka Spur (95 mW m -2 ) was also likely produced by the topographic effect because the measurement position is close to the steep western slope of the Umitaka Spur. Variable surface heat flow values (81-105 mW m -2 ) were obtained at and around the cold seep sites in the summit area on the Umitaka Spur. Local advective heat transport owing to local pore fluid flow in the sediment associated with cold seep activity might result in such variation of surface heat flow values.

Estimation of the regional background heat flow in the Joetsu Basin
Surface heat flow data obtained in this study do not represent the regional background heat flow from deep in the basin owing to the topographic effect (Fig. 8). Lachenbruch (1968) presented a 2D analytical solution for calculating surface heat flow disturbed by topographic change that can be used to correct for the topographic effect.
However, the correction is limited to areas with simple topography and the accuracy is insufficient for calculating surface heat flow disturbed by complicated topography (e.g., Ganguly et al. 2000;He et al. 2007).
To investigate the regional background heat flow in the Joetsu Basin, we created a 2D finite element geothermal structure model that involves the topographic change and subsurface structure of the METI 2D seismic survey line in the basin (Figs. 2 and 8). The model was used to calculate the distribution of topographically disturbed surface heat flow at the seafloor by assigning various heat flow values to the bottom boundary of the model. COMSOL Mul-tiphysics® software (COMSOL AB, Stockholm, Sweden) was used to calculate the geothermal structure and surface heat flow. The regional background heat flow in the Joetsu Basin was estimated by comparing the calculated surface heat flow distribution and measured surface heat flow values. The surface heat flow data obtained at station HF1 were not used because of the possible sedimentation effect, as discussed previously. The surface heat flow data likely undisturbed by BWT fluctuations obtained by Machiyama et al. (2009) were also not used because the topography of the measurement positions differed from those in the northern to northeastern part of the Umitaka Spur where the METI 2D seismic survey line is located (Fig. 2). Figure 9 shows the 2D model used in the present study to calculate the geothermal structure and surface heat flow based on the topography and structure of the METI 2D seismic profile (Fig. 8b). A seawater P-wave velocity of 1474 m s -1 is assumed to calculate the water depth (Japan Oil, Gas and Metals National Corporation 2005). A timedepth conversion curve determined by vertical seismic profiling conducted at wells A and B (Japan Oil, Gas and Metals National Corporation 2005) was used to convert the two-way travel times to a depth below the seafloor for each formation boundary. No information is available regarding the bottom boundary of the Green Tuff along the seismic profile. We therefore assumed that the Green Tuff occurs below the Nanatani Formation and the bottom boundary of the model was set to be sufficiently deep (Fig. 9b). The finite elements for the geothermal structure modeling are shown in Fig. 9c.

Structural model and boundary conditions
For the model boundary conditions, the seafloor was fixed as the top boundary with a constant temperature of 0 °C, and both sides of the model were assumed to be adiabatic. Various heat flow values were set to the bottom boundary of the model as the regional background heat flow.

Physical properties model
The thermal conductivity of sediment to a depth of 34 m below the seafloor in the Joetsu Basin was measured from sediment core samples recovered during the MD179 cruise (Goto et al. 2017). However, thermal conductivity data are not available below 34 m. The thermal conductivity of sediment for the geothermal structure modeling was therefore calculated using existing porosity data and an empirical relationship between the porosity and thermal conductivity of sediment. Figure 10 shows the vertical porosity distribution measured from samples recovered during the MD179 cruise (Goto et al. 2017) and at well C (Toyama 1985). The porosity data by Goto et al. (2017) were obtained to a depth of 36 m below the seafloor. Porosity data were obtained at well C at a depth interval between 572 and 3081 m below the seafloor (Fig. 10). Thus, porosity data are presently unavailable for depths 36-572 m below the seafloor. The dominant lithology of the Haizume Formation and Nishiyama Formation is mud. We assumed that the porosity of the formations decreases exponentially with depth and found that the vertical porosity distribution can be modeled by where z 1 is the depth from the seafloor and (z 1 ) is the porosity of sediment at depth z 1 . The modeled porosity distribution for the Haizume Formation and Nishiyama Formation is plotted as a thick solid line in Fig. 10. For the Shiiya, Lower Teradomari, and Nanatani formations, which are composed of mud, sand, and tuff, we applied a porosity model by fitting the porosity data of the formations (dashed curve in Fig. 10) to where z 2 is the depth from the top of the Shiiya Formation and (z 2 ) is the sediment porosity at depth z 2 . The Upper Teradomari Formation has a dominant lithology of mud (Muramoto et al. 2007) and occurs between the Shiiya and Lower Teradomari formations. We used the porosity model (Eq. (4)) to estimate the porosity of the Upper Teradomari Formation.

Porosity model
(3) z 1 = 0.718e −1.07×10 −3 z 1 (4) z 2 = 0.237e −2.56×10 −5 z 2 Fig. 9 Two-dimensional structural model for geothermal structure and surface heat flow calculations. a Magnified view of the structural model. Boundaries between the Haizume Formation and Nishiyama Formation and between the Lower Teradomari Formation and Nanatani Formation are indicated by thin dashed lines. b Full view of the structural model. c Finite elements for geothermal structure modeling Fig. 10 Vertical porosity distribution measured from samples recovered during the MD179 cruise (Goto et al. 2017) and at well C (Toyama 1985). Dotted horizontal lines show the formation boundaries. Solid and dashed curves represent the porosity models for the calculation of thermal conductivity in the formations 1 3 Porosity data are also not available for the Green Tuff. We therefore assumed that the Green Tuff has uniform porosity and used an average porosity value (11.1%) of rock samples from the Green Tuff regions in eastern Japan (Deguchi et al. 1996).

Thermal conductivity model
The thermal conductivity of the subsurface material for geothermal structure modeling was calculated from the porosity model (Fig. 10) using an empirical relationship between porosity and sediment thermal conductivity following the geometric mean model, as proposed by Woodside and Messmer (1961). The geometric mean model approximates the relationship between the porosity and thermal conductivity of sediment using the following formula: where K s is the thermal conductivity of sediment, is the sediment porosity, K w is the thermal conductivity of seawater that fills the pore spaces in sediment, and K g is the thermal conductivity of the sediment grains. Goto et al. (2017) estimated K g to be 2.76 W m -1 K -1 by assuming a thermal conductivity of seawater (K w ) of 0.596 W m -1 K -1 (Kaye and Laby 1995)) and fitting Eq. (5) to the porosity and thermal (5) K s = K w K 1− g conductivity data for mud-dominant sediment in the Joetsu Basin. This K g value was therefore also used to calculate the thermal conductivity in the Haizume, Nishiyama, and Upper Teradomari formations in the present study.
The Shiiya, Lower Teradomari, and Nanatani formations have lithologies consisting of an aggregate of mud, sand, and tuff. Thermal conductivity data are not available for sediment grains of this lithology. The K g of the formations was therefore assumed to be an arithmetic mean of the K g values of mud-dominant sediment (2.76 W m -1 K -1 ; Goto et al. 2017) and sandy sediment (3.54 W m -1 K -1 ), which was estimated for sandy sediment in the eastern flank of the Juan de Fuca Ridge by Goto and Matsubayashi (2009).
The K g of the Green Tuff was assumed to be 3.59 W m -1 K -1 , which is an average value reported for rock samples obtained in the Green Tuff regions in eastern Japan (Deguchi et al. 1996).
We applied the Sekiguchi model (Sekiguchi 1984) to account for the temperature dependence of K g for any lithology. The thermal conductivity of the pore fluid in sediment was assumed to be 0.596 W m -1 K -1 .

Regional background heat flow in the Joetsu Basin
We calculated surface heat flow from the geothermal structure model by setting various constant heat flow values at Fig. 11 Surface heat flow distribution and geothermal structure along the METI 2D seismic profile. a Surface heat flow values measured during the MD179 cruise and surface heat flow distributions calculated from the geothermal structure models with three regional background heat flow values, which well describe the measured values.
Station names for surface heat flow measurements are indicated. b Geothermal structure for an assumed regional background heat flow value of 85 mW m -2 . Isotherms are indicated as thin solid lines. Thick solid and thin dashed lines indicate the formation boundaries shown in Fig. 9 the bottom boundary of the model and visually searching the range of regional background heat flow values that can explain those obtained during the MD179 cruise. We found that the measured surface heat flow distribution can be explained by a regional background heat flow value of 85 ± 6 mW m -2 (Fig. 11a). Figure 11b shows the geothermal structure calculated using a background heat flow of 85 mW m -2 . The geothermal structure near the seafloor is controlled by the seafloor topography; for example, a 50 °C isotherm is sub-parallel to the seafloor. The topographic effects on the geothermal structure decrease with increasing depth below the seafloor. The thermal structure becomes increasingly influenced by the geological structure with depth. The 100 and 150 °C isotherm depths increase from southeast to northwest. In contrast, the 200 °C isotherm depth increases from northwest to southeast. These isotherms are likely controlled by the thickness of the sedimentary rock above the Green Tuff and the top depth of the Green Tuff, which has a higher thermal conductivity than the overlying sedimentary rock.
We evaluated the effects of the assumed sediment thermal conductivity model on the calculation of the geothermal structure model and surface heat flow using K g values ± 10% of the originally assumed values. The evaluation indicates that the geothermal structure model indeed depends on K g : applying 10% lower K g resulted in a higher geothermal structure model than that calculated using the originally assumed values. The temperature difference between these geothermal structure models increases with depth, reaching ~ 13.1 K on the top of the Green Tuff on the northwestern side of the model. The opposite results were obtained for K g values 10% higher than the originally assumed values. In this case, the calculated geothermal structure was lower than the originally calculated structure, and the temperature difference reached ~ 11.6 K on the top of the Green Tuff on the northwestern side of the model. However, surface heat flow values calculated from the geothermal structure models using K g values ± 10% of the originally assumed values were nearly the same as those calculated from the originally calculated geothermal structure model. This indicates that the thermal conductivity models have less of an effect on the estimation of regional background heat flow.
We calculated the regional background heat flow, geothermal structure, and surface heat flow expected from the geothermal structure for the Joetsu Basin. These provide geothermal constraints on other geophysical, geochemical, biological, and geological factors in the basin. For example, basin modeling is required to infer petroleum systems (e.g., Hantschel and Kauerauf 2009). The regional background heat flow provides a constraint on the present-day heat flow for constructing historical heat flow models of the basin to infer the generation, migration, and accumulation of oil and gas. The geothermal structure can be used to investigate gas hydrate stability (e.g., Dillon and Max 2000), subseafloor microbial ecosystems (Heuer et al. 2020), and the present-day generation depths of oil and gas (Hunt 1996). In areas with locally upward pore fluid flow in sediment, higher surface heat flow can be expected to be measured. However, higher relative surface heat flow can be measured in areas with concave topographic changes (Lachenbruch 1968;Blackwell et al. 1980), even if upward pore fluid flow does not occur in the sediment. Surface heat flow values expected from a geothermal structure model that considers topographic changes provide a baseline surface heat flow for investigating advective heat transport owing to local pore fluid flow in sediment.

Surface heat flow measurement in areas with significant BWT fluctuations
In the present study, the depth at which the effect of BWT fluctuations on the sediment temperature decays was evaluated using long-term BWT fluctuation data. Surface heat flow values undisturbed by the BWT fluctuations were successfully determined from sediment temperatures measured at depths deeper than the critical depth using a 15 m long geothermal probe.
There are two requirements to accurately measure surface heat flow using a long geothermal probe in marine areas with significant BWT fluctuations. One is that the penetration of a long geothermal probe depends on the lithology and consolidation state of the sediment. For example, increasing the volume fraction of sand in sediment impedes the probe penetration, and the geothermal probe cannot penetrate consolidated sediment. Surface heat flow measurements using a long geothermal probe are possible in areas covered by soft sediment, such as the mud-dominant sediment in the Joetsu Basin. The other requirement is that information regarding BWT fluctuations and sediment thermal diffusivity in the survey area are required to evaluate the depth at which the sediment temperature is thermally affected by BWT fluctuations. This depth depends on the amplitude and frequency of the BWT fluctuations and sediment thermal diffusivity (Carslaw and Jaeger 1959;Goto et al. 2005). The information can also be used for the special design of instruments, such as an optimized geothermal probe length, to measure surface heat flow undisturbed by BWT fluctuations.
Effect of thermal disturbance owing to topographic changes using a 2D geothermal structure model Surface heat flow values obtained along a 2D seismic survey line in the Joetsu Basin showed thermal disturbances owing to topographic changes (Fig. 11). We estimated the regional background heat flow in the basin by comparing the surface heat flow distributions calculated from a 2D geothermal structure model using various constant heat flow values at the lower boundary and measured heat flow values at the seafloor (Fig. 11).
The effect of topography on surface heat flow depends on the complexity of the topographic change. He et al. (2014) determined BSR-derived heat flow from three-dimensional (3D) seismic data obtained in the mid-slope region of northern Cascadia subducting margin. They used 2D and 3D geothermal models to correct for the topographic effect on BSR-derived heat flow and revealed that topographic disturbances owing to complicated topographic changes cannot be corrected using a 2D geothermal model and a 3D geothermal model is required. They also indicated that gentle topographic disturbances have a small effect on surface heat flow and the difference between topographic disturbances calculated by 2D and 3D geothermal models is small.
The METI 2D seismic survey line used for our 2D geothermal modeling was set nearly perpendicular to the long axis of the Joetsu Knoll (Fig. 2). The topography of the Joetsu Knoll along the seismic profile is relatively simple (Fig. 8). We therefore consider this situation to be reasonable for use of a 2D geothermal structure model to investigate regional background heat flow. Surface heat flow values measured at stations HF4-HF7 on and around the Joetsu Knoll can be explained using geothermal structure models with a regional background heat flow of 85 ± 6 mW m -2 (Fig. 1112).
The METI 2D seismic survey line passes through terrace-like topography on the northern to northeastern slope of the Umitaka Spur (Fig. 2), which is characterized by gentle topography. The surface heat flow values measured at stations HF2 and HF3 can also be explained using a 2D geothermal model with a regional background heat flow of 85 ± 6 mW m -2 (Fig. 11), which indicates the small effect of topographic changes on the surface heat flow values at those stations.

Conclusion
To obtain reliable surface heat flow data in the Joetsu Basin in the eastern margin of the Japan Sea where BWT fluctuations have been suggested to disturb sediment temperatures near the seafloor, we measured sediment temperatures at seven stations in the basin by penetrating a 15 m long geothermal probe to a depth where the effect of BWT fluctuations is sufficiently negligible. Our results of a long-term BWT monitoring study indicate that BWT fluctuations can disturb sediment temperatures at depths less than 2 m below the seafloor. Surface heat flow values of 81-88 mW m -2 were determined from undisturbed sediment temperature data measured between 2.3 and 16.3 m below the seafloor.
A comparison of the obtained surface heat flow data and topography from an existing 2D seismic profile reveals that surface heat flow values at six stations show thermal disturbances owing to topographic changes. We constructed a 2D geothermal structure model to reflect the topography and geologic structure along the 2D seismic profile and estimated a regional background heat flow of 85 ± 6 mW m -2 by comparing the measured surface heat flow values with the surface heat flow distributions calculated from the 2D geothermal structures.
Two requirements are proposed to accurately measure surface heat flow undisturbed by BWT fluctuations in marine areas using a long geothermal probe: (i) information regarding the lithology and consolidation state of the seafloor sediment for probe penetration; and (ii) information regarding the BWT fluctuation and thermal diffusivity of sediments to accurately evaluate the depth to which BWT fluctuations affect the sediment temperature. If these requirements are not satisfied, other techniques must be applied to obtain heat flow data undisturbed by BWT fluctuation, such as heat flow measurements in a deep borehole, long-term monitoring of bottom-water and sediment temperatures to obtain sediment temperatures for which the effects of BWT fluctuations are corrected, and heat flow estimations from BSR depths if a where T(z) is the temperature at depth z, T 0 is the temperature at z = 0, q is the surface heat flow, and R(z) is the thermal resistance expressed by In the present study, K(z) is approximated as where a, b, and m are constants estimated by fitting the sediment thermal conductivity data obtained in the Joetsu Basin (Fig. 12). Using Eq. (A3), the thermal resistance R(z) becomes Uncertainties of the sensor depth and modeled thermal conductivity of the sediment result in thermal resistance uncertainties. In the present study, the thermal resistance uncertainties are calculated by error propagation (Taylor 1997): where σ R is the standard deviation of the thermal resistance R, σ z is the standard deviation of the sensor depth, and σ a , σ b , and σ m are the standard deviations of a, b, and m in Eq. (A3), respectively. In the present study, the sensor depth uncertainties (see surface heat flow measurement section in the main text) estimated from the standard deviation of the BWT fluctuations were used as the standard deviation of the sensor depth σ z .
To estimate the surface heat flow from sediment temperatures and thermal resistance using Eqs. (A1) and (A4), we applied the weighted orthogonal distance regression that accounts for the sediment temperature and thermal resistance uncertainties using Igor Pro software (WaveMetrics, Inc., Lake Oswego, OR, USA).
Data availability Because this study was conducted as part of a national project of methane hydrate research in the Japan Sea, we do not have permission to share the data obtained in this project. We also do not have permission to share the METI two-dimensional seismic data and the wells A and B data obtained by other national research projects in Japan.

Conflicts of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.