A Long-Term Geothermal Observatory Across Subseafloor Gas Hydrates, IODP Hole U1364A, Cascadia Accretionary Prism

We report 4 years of temperature profiles collected from May 2014 to May 2018 in Integrated Ocean Drilling Program Hole U1364A in the frontal accretionary prism of the Cascadia subduction zone. The temperature data extend to depths of nearly 300 m below seafloor (mbsf), spanning the gas hydrate stability zone at the location and a clear bottom-simulating reflector (BSR) at ∼230 mbsf. When the hole was drilled in 2010, a pressure-monitoring Advanced CORK (ACORK) observatory was installed, sealed at the bottom by a bridge plug and cement below 302 mbsf. In May 2014, a temperature profile was collected by lowering a probe down the hole from the ROV ROPOS. From July 2016 through May 2018, temperature data were collected during a nearly two-year deployment of a 24-thermistor cable installed to 268 m below seafloor (mbsf). The cable and a seismic-tilt instrument package also deployed in 2016 were connected to the Ocean Networks Canada (ONC) NEPTUNE cabled observatory in June of 2017, after which the thermistor temperatures were logged by Ocean Networks Canada at one-minute intervals until failure of the main ethernet switch in the integrated seafloor control unit in May 2018. The thermistor array had been designed with concentrated vertical spacing around the bottom-simulating reflector and two pressure-monitoring screens at 203 and 244 mbsf, with wider thermistor spacing elsewhere to document the geothermal state up to seafloor. The 4 years of data show a generally linear temperature gradient of 0.055°C/m consistent with a heat flux of 61–64 mW/m2. The data show no indications of thermal transients. A slight departure from a linear gradient provides an approximate limit of ∼10−10 m/s for any possible slow upward advection of pore fluids. In-situ temperatures are ∼15.8°C at the BSR position, consistent with methane hydrate stability at that depth and pressure.


INTRODUCTION
The stability of methane hydrates in continental margin sediments has long been known to be strongly dependent on temperature and pressure (e.g., Kvenvolden and Barnard, 1983;Hyndman et al., 1992), such that seismically determined depths to bottom-simulating reflectors (BSRs) that often mark the base of gas hydrate stability (BGHS) have been used by many authors to estimate heat flux through the sedimentary section above the BSR (e.g., Yamano et al., 1982;Davis et al., 1990;Hyndman et al., 1992). However, there have been only a few direct measurements of temperature within the gas hydrate stability zone (GHSZ) made in scientific ocean drilling, and even fewer long-term time-series measurements of temperature across the GHSZ. The best previous example of the latter was at Ocean Drilling Program (ODP) Site 892 at Hydrate Ridge in the southern Cascadia accretionary prism, where thermal results were enigmatic with respect to inferred depths of methane hydrate stability (Davis et al., 1995). This paper describes 4 years of temperature time-series data in the northern Cascadia accretionary prism, at Integrated Ocean Drilling Program (IODP) Site U1364, spanning a clear BSR and the GHSZ. These data are consistent with and resolve more clearly the previously interpreted temperature gradient and temperature at the BGHS at the site. They also improve the estimate of heat flux at depth and provide an upper limiting constraint on the rate of any vertical pore-fluid expulsion.
As shown in Figure 1, Site U1364 lies ∼20 km landward of the toe of the Cascadia subduction zone accretionary prism, where much of the thick section of turbidite and hemipelagic sediments deposited on the eastern flank of the Juan de Fuca Ridge is scraped off the underthrusting oceanic crust (Davis and Hyndman, 1989;Hyndman et al., 1992;Westbrook et al., 1994). Convergence of the Juan de Fuca oceanic plate relative to the North American continental plate occurs in a direction roughly normal to the continental margin at a rate of about 42 mm/yr (DeMets et al., 1990). A topographic trench at this subduction zone is absent as a consequence of the extremely high rate of turbidite sediment supply from the continent during the Pleistocene (Davis and Hyndman, 1989). At the accretionary prism toe or deformation front, where the seawardmost thrust fault of the accretionary complex intersects the seafloor, the sediments that bury the eastern Juan de Fuca Ridge flank are approximately 2.7 km thick. At Site U1364 the accreted sedimentary section is nearly doubled to a thickness of approximately 5 km (Yuan et al., 1994). With such tectonic thickening and compaction, pore fluids are expelled, and gas -primarily methane -is transported upward to contribute to the formation of gas hydrates in the upper few hundred meters of the sediment section (Davis et al., 1990;Hyndman and Davis, 1992;Haacke et al., 2007;Riedel et al., 2010a).
Site U1364 lies at a position where the fluid expulsion rate, estimated on the basis of the rates of compaction and vertical growth of the prism, reaches a cross-margin maximum, and where a clearly developed BSR marks the BGHS (Davis et al., 1990;Hyndman and Davis, 1982). That rationale had previously been used in selecting the locations of adjacent ODP Site 889 and IODP Site U1327; the three sites essentially share a common location within a diameter of 700 m ( Figure 2), with each being at least 1 km away from any mapped locations of seafloor cold vents (Scherwath et al., 2019). Cores and logs from Sites 889, 890, and U1327 Riedel et al., 2010a) have documented the nature of the incoming Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 568566 undeformed sediments, the compaction history during accretion, the details of the lithology, and the distribution and composition of gas hydrates across the area. In particular, they show that 1) pore fluids are significantly freshened below the upper unit of undeformed slope sediments, which should affect estimates of the depth limit of the GHSZ; and 2) there are geochemical and carbon isotope signatures for primarily microbial gas generation at this location, as summarized by Riedel et al. (2010a). At the three previous sites, ten subseafloor sediment temperature measurements down to ∼210 mbsf constrained an average temperature gradient of 0.055°C/m, with an estimated temperature of 15-16°C at the BGHS estimated at ∼230 mbsf in the area (Exp 311 Scientists et al., 2006a;Riedel et al., 2010a).

METHODS: BOREHOLE OBSERVATORY AND TEMPERATURE SENSOR CONFIGURATION
During IODP Expedition 328 in 2010, Hole U1364A was drilled without coring or any logging to a total depth of 336 mbsf through roughly 90 m of gently deformed slope-basin deposits and underlying sediments of the accretionary prism that are folded and faulted on a scale too small to be resolved in seismic reflection profiles ( Figure 3). An "Advanced CORK" (ACORK) borehole observatory was installed immediately after drilling (Davis et al., 2012), configured as shown in Figure 3A projected onto trench-normal seismic line PGC 89-08. The ACORK was constructed around solid 10.75-inch casing, left open at the wellhead above seafloor but sealed at the bottom of the casing with a bridge plug backed with cement. This left the interior open to a depth of 302 mbsf for instruments requiring thermal or mechanical-but not direct-contact with the formation. Formation fluid pressures are transmitted to sensors at the wellhead via 2.03-m-long circumferential sand-packed filter screens and 3-mm-diameter (1/8-inch i.d.) stainless steel hydraulic tubing mounted on the outside of the casing. Screens (numbered S1 to S4 from bottom to top) are positioned at depths of 304 mbsf, 244 mbsf, 203 mbsf, and 156 mbsf. The middle two pressure monitoring screens were positioned 14 m below and 27 m above the base of the gas hydrate stability zone at 230 mbsf to observe the effects of free gas and gas hydrate in the sediment matrix, and diffusive signals originating at the hydrate-gas boundary. The lowermost and uppermost screens were placed 74 m below and above the boundary, a distance that was anticipated to be sufficient to avoid hydrologic complications originating at the boundary. Initial results (Davis et al., 2012) indicated that collapse of the formation around the installation provided good hydraulic isolation of the three deepest screens, but the seafloor reference pressure gauge had failed. Therefore, the input to the gauge originally assigned to the uppermost screen was switched from formation to the seafloor to provide a continuing reference pressure signal. In July of 2016, we installed a combined seismometer/ tiltmeter/thermistor array in the central bore of the ACORK, using the ROV Jason from R/V Sikuliaq. An earlier depth check conducted in May of 2014 with a Seabird conductivitytemperature-depth (CTD) profiler on an ROV-mounted winch had indicated the likely presence of fill below 300 mbsf. To minimize risk on deployment, the 2016 sensor string was designed to terminate above that section, somewhat above the deepest pressure-monitoring screen. The seismometer/tiltmeter package was mechanically clamped to the casing at 280 mbsf and cabled to a master control unit that landed in the ACORK wellhead (McGuire et al., 2018). During deployment, a 24thermistor cable was married to the seismometer/tiltmeter string, with the deepest thermistor positioned at 268 mbsf, roughly midway between the two deepest screens. The thermistor cable was a Concerto T24 unit manufactured by Richard Brancker Research Ltd., with nominal resolution of 0.001°C in the expected range of 2-20°C. It had been designed with variable spacing between sensors of 4-18 m, with closely spaced sections around the upper three pressure monitoring screens and BSR and wider spacing elsewhere ( Figure 3B). The temperature sensors are numbered T1 to T24 from bottom to top.
The output from the temperature data logger was fed to the master sensor string control unit, which was also designed to Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 568566 allow direct connection for power and data streaming to the Ocean Networks Canada (ONC) cable at its Clayoquot Slope node. That connection could not be completed in 2016 but was successfully made in June of 2017, after which temperature data were logged continuously at one-minute intervals with the exception of a few short ONC outages. During the 2016 deployment cruise, about 1 day of initial temperature data had been logged to a shipboard computer while we remained connected to the installation via Jason. Unfortunately, in May of 2018 the main internet switch in the master control unit failed, ending the streaming of all data including temperature. The entire 2016 inside-casing instrument string was pulled from the hole in Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 568566 4 late summer 2019, also using Jason from R/V Sikuliaq, for replacement of the internet switch, refurbishment, and redeployment hopefully in the near future. Monitoring of pressure data via an independent ONC cable connection has not been affected.

RESULTS
Hole U1364A Temperature Time-Series Data Figure 4 shows the eleven-month time-series of thermistor cable data from the June 2017 ONC connection to the May 2018 failure of the internet switch in the main control unit for the sensor string. For most of the sensors, the long-term records are quite stable, but four or five sensors displayed indications of sudden or gradual degradation, possibly due to seawater penetration into the cable. For four of those sensors, only the more stable early sections of data are shown on Figure 4 and used in subsequent analyses. In particular, readings from T14 (154 mbsf) and T10 (200 mbsf) showed abrupt deviations and subsequent unstable readings starting in late September 2017 and mid-March 2018, respectively. In addition, readings from T19 (76 mbsf) and T5 (232 mbsf) showed upward drifts with steadily increasing noise levels (to >1°C) beginning in August and October 2017, respectively. The record from T2 (256 mbsf) also shows a slow upward drift but without increasing noise levels, eventually reaching higher apparent temperatures than the deepest sensor T1 (268 mbsf). This is probably due to slow failure of T2, but we cannot definitively rule out some other cause such as formation warming at that depth or slow convective overturn deep in the hole. Convection of the fluids within casing is not unexpected for a hole of this diameter and temperature gradient, but it does not normally disturb the average gradient (e.g., Diment, 1967;Sammel, 1968). Therefore, the T2 readings are still shown on Figure 4, but they were not used in analyses of temperature gradient or heat flux in sections below. In addition, the record from T24, at a position within the ACORK but above seafloor, clearly shows tidal variations in bottom-water temperatures, and the T23 record from just below seafloor shows more attenuated tidal variations that are barely discernible on The best-fit linear gradient from our data of 0.055°C/m is consistent with the value reported for the average of ODP/IODP sediment probe readings at nearby Sites 889, 890, and U1327 (Expedition 311 Scientists et al., 2006b). Despite this consistency, Figure 5 shows that several of those ODP/IODP probe readings seem to plot below our profiles, especially data from Sites 889/ 890. This may be due to a combination of temporal variability of bottom water temperatures as discussed below in Section 4.3 and uncertainties introduced by the normal ODP/IODP practice of correcting probe readings so that temperatures within the drillpipe measured during a stop at mudline depths match known bottom-water temperature. That correction is based on the assumption that the drillpipe acts as a perfect heat exchanger, cooling the surface seawater used to pump the probe down the pipe to the bottom-water value. The corrections for this effect applied to the two deepest readings from Site 889 shown on Figure 5 were substantial, −1.3°C (Shipboard Scientific Party, 1994), so there may be larger uncertainties with those data compared to our 2014-2018 data. The 2014-2018 profiles shown in Figure 5 show no clear indications of transients at the level of the BSR or pressuremonitoring screens during the four-year period of collection of temperature data in Hole U1364A. The profiles do hint at a small degree of convex-upward curvature that might arise from diffuse upward pore-fluid migration associated with the gas hydrate stability zone as modeled by Hyndman and Davis (1982) and Riedel et al. (2010a). Estimating heat flux at depth and properly assessing an upper limit on possible fluid advection consistent with the temperature profiles first requires careful assessment of the in-situ thermal conductivity.

Thermal Conductivity and Heat Flux Estimates
Unfortunately, there are much greater uncertainties in the in situ thermal conductivity than the temperature profile. As reported by Shipboard and Scientific Party (1994), Expedition 311 Scientists et al. (2006b), and Riedel et al. (2010b), thermal conductivity measurements on ODP and IODP cores from Sites 889 and U1327 display considerable scatter with no clear depth trend. Much of the scatter might be due to variable disturbances in the coring and recovery processes, including de-gassing effects. In addition, Riedel et al. (2010b) noted that Site 889 values are about 10% higher on average than Site U1327 values. The same needleprobe method was used (Von Herzen and Maxwell, 1959), but with different instruments; this is unlikely the cause of the difference in conductivities, but an uncertainty of at least 10% cannot be eliminated (Riedel et al., 2010b). Therefore, Expedition 311 Scientists et al. (2006b) and Riedel et al. (2010b) used an average value of 1.1 W/m-K ± 10% to estimate a heat flux of 61-62 mW/m 2 at the sites. Applying that same conductivity value to our nearly linear gradient of 0.055°C/m yields a similar estimate of 60.5 mW/m 2 conductive heat flux through the section.
The Expedition 311 Scientists et al. (2006b) also noted that the conductivities for the more competent cores approached the "regional trend" of conductivity vs. depth developed by Davis et al. (1990) based on a seismic velocity to porosity to conductivity transform. Where thermal conductivity varies in a known way with depth, Bullard (1939) first suggested a transform to thermal resistance, the integral to a given depth of the inverse of the conductivity, with the slope of a plot of temperature vs thermal resistance providing a determination of conductive heat flux. The Davis et al. (1990) computation suggested a variation of conductivity that could be fit by a second order polynomial in depth with a seafloor value of 1.07 W/m-K. In an independent and potentially more reliable approach specific to the site, we also applied the relationship between thermal conductivity and porosity developed by Goto and Matsubayashi (2009) for Cascadia Basin sediments to the logging-while-drilling (LWD) porosity profile collected in Hole U1327A. The LWD data showed porosities decreasing with depth down to ∼120 mbsf, a higher porosity zone at ∼120-140 mbsf, and more uniform average values in the deeper section other than higher values just below the BSR (Expedition 311 Scientists et al., 2006b). This approach allowed computation of a detailed profile of the variation of in situ thermal conductivity with depth, and an integration to estimate thermal resistance at each thermistor depth. As shown in Figure 6, plots of thermistor cable temperatures vs. thermal resistance are quite linear for both models. The slopes yield conductive heat flux values of 63 mW/m 2 using the Davis et al. (1990) model and 64 mW/m 2 using the Goto and Maysubayashi (2009) model.

Constraints on Advection Rates
The linearity of the temperature-thermal resistance plots in Figure 6 suggests that an increase of conductivity with depth might explain at least some of the slightly convex-upward nature of the temperature-depth plots in Figure 5. Given the uncertainties in the depth variation of the in-situ conductivity profile, it is debatable whether it is valid to interpret an upward advection rate from the slight apparent curvature in the temperatures vs depth. Nevertheless, we applied the method of Bredehoeft and Papadopulos (1965) to define an upper limit to a possible advection rate consistent with the temperature data. That analysis suggests a value of 0.1 for the Bredehoeft and Papadopulos (1965) "β" parameter over the 268 m subseafloor depth range of the May 2018 thermistor cable data (Figure 7). β is defined as vρcL/k, where v is vertical advection rate, ρ and c are density and specific heat of the pore fluid, respectively, L is the vertical range of temperature measurements (268 m), and k is average thermal conductivity in that range. Using 1.1 W/m-K for average conductivity, 1,030 kg/m 3 for density of somewhat freshened pore fluids, and 3850 J/kg-K for specific heat of the pore fluids yields an estimate of upward vertical advection rate on the order of ∼1 × 10 −10 m/s. This is remarkably similar in order of magnitude to the maximum expulsion rate of nearly 3 mm/yr estimated by Hyndman and Davis (1982) at the general position of Hole U1364A on the northern Cascadia margin and background methane flux rates estimated by Riedel et al. (2006a) at Site U1328. We consider our estimate an upper limit for reasons described above and because the Bredehoeft and Papadopulos (1965) method is especially sensitive to the value of the deepest temperature reading, which in this case falls below the generally linear trends shown in Figures 5 and 6.

Relevance to Regional Heat Flux Measurements
Several hundred km farther south along the Cascadia Margin, Tréhu (2006) conducted an extensive program of downhole temperature probe measurements during ODP Leg 204 drilling at southern Hydrate Ridge offshore Oregon. That study produced comparable results to ours: a slightly lower average heat flux value of ∼55 mW/m 2 and nearly linear temperature-depth profiles generally consistent with predicted temperatures at the base of the GHSZ as estimated from BSR depths there. As reported in Section 4.1, the 61-64 mW/m 2 range of heat flux values we obtain from the thermistor cable readings in Hole U1364A is very close to the values obtained in nearby Sites 889, FIGURE 6 | Thermistor cable temperatures in Hole U1364A plotted against integrated thermal resistance, using the models of Davis et al. (1990) and Goto and Matsubayashi (2009) for the depth variation of thermal conductivity at the site. Note that the temperature scales for the two plots are offset by 2°C to avoid direct overlap.
FIGURE 7 | Non-dimensionalized May 2018 thermistor cable readings plotted against the Bredehoeft and Papadopulos (1965) f(β,z/L) function along with characteristic type curves values for values of β associated with upward migration of pore-fluids. As detailed in the text, the consistency of the temperature data with the type curve for β 0.1 suggests an upward porefluid migration rate no greater than ∼10 −10 m/s. Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 568566 890, U1327, and U1328. It is also close to the average value of 58 mW/m 2 obtained with a 4-m long probe during a regional heat flux survey in 2000 with trench-normal and trench-parallel lines crossing over Site U1328 . However, earlier trench normal heat flow surveys with a 2.5-m long probe in the same area yielded significantly higher values (Davis et al., 1990). Riedel et al. (2006b) suggested that this discrepancy might be attributed to time-varying perturbations to shallow sediment temperatures from bottom-water temperature changes at time scales of a few months to years. This interpretation is supported by the combination of our significantly deeper heat flux determination in Hole U1364A, and a ten-year bottom-water temperature time-series obtained with a nearby bottom pressure recorder (BPR) ∼3 km away near the ONC Clayoquot Slope node (Figure 2). Daily running averages of the BPR temperatures from 2015-2020 are plotted in Figure 8 and clearly show variations in bottom-water temperature of up to ±0.15°C over time scales of months to nearly 2 years. The analysis of Davis et al. (2003) shows that downward propagation of bottom-water temperature variations of those magnitudes and time scales into the sediments will produce significant temporal perturbations to temperature gradients in the uppermost 2-3 m of sedimentsthe depth range penetrated by shallow heat flow probes. If similar variations in bottom-water temperatures were occurring at the times of the early heat flux surveys, they could indeed explain the differences between the Davis et al. (1990), Riedel et al. (2006a), and Riedel et al. 2006b) regional heat flux trends. This argues for caution in using short-probe heat flux determinations on the upper slopes of accretionary prisms as thermal constraints on models of deep-seated tectonic processes and state at subduction zones.
Deeper borehole data like the thermistor cable data from Hole U1364A -and downhole probe data like those of Tréhu (2006) are out of range of the effects of such bottom-water temperature perturbations, and in that sense provide much more reliable measurements of seafloor heat flux and deep thermal structure.
While our data and those of Tréhu (2006) represent just two well documented studies, it is notable that the downhole temperature profiles are consistent with predicted methane hydrate stability at both locations. This further validates the technique of estimating seafloor heat flux from BSR depths, as long as BSR depths are well constrained by seismic data and thermal conductivity profiles can be accurately determined.

CONCLUSION
Four years of discrete and continuous temperature logging in Hole U1364A on the Vancouver Margin accretionary prism has provided important constraints on the geothermal state and an upper limit on possible pore-fluid flow in the prism associated with the BSR and BGHS penetrated at ∼230 mbsf at this location. The temperature profiles have remained quite constant over those 4 years, defining a nearly linear temperature gradient of 0.055°C/ m with no indications of in-situ transients during the period. This gradient is consistent with a conductive heat flux of 61-64 mW/ m 2 depending on the model adopted for the less constrained variation of thermal conductivity with depth at the site. The generally linear temperature profiles could also be consistent with a slight convex-upward curvature associated with upward porefluid migration at low rates no greater than 10 −10 m/s, the same order as previously modeled by Hyndman and Davis (1982) to bring methane from depth to the gas hydrate stability zone. Finally, the temperature data show in-situ values at the BSR depth of ∼15.8°C, consistent with predicted methane hydrate stability at that depth and pressure.

DATA AVAILABILITY STATEMENT
The thermistor cable temperature data used in this report are available from the Ocean Networks Canada data portal as specified in Ocean Networks Canada Society (2020). (Note that the archived thermistor cable data include periodic spikes during the first 6 months of 2017 recording that were caused by regular ground-fault checks of the WHOI borehole instruments; those spikes have been filtered out from any data used in this report.) The 2014 CTD data from Hole U1364A are available on Scholars Portal Dataverse as specified in Davis et al. (2020).

AUTHOR CONTRIBUTIONS
KB obtained NSF funding for the thermistor cable installed in Hole U1364A, participated in the design and installation of the combined sensor string installed in 2016, and led interpretation of the thermistor cable data and writing of this manuscript. ED was chief scientist of the IODP expedition that installed the ACORK in Hole U1364A in 2010, led interpretation of pressure data from the ACORK and the 2014 effort to conduct a depth check in the hole with CTD, and collaborated on the design of the thermistor cable and interpretation of the geothermal data. MH participated in the ACORK installation expedition and interpretation of 8 | Daily running averages of the internal temperature sensor for BPR NC89 deployed near the ONC Clayoquot Slope node, ∼3 km southeast of Hole U1364A. Also shown is the 11-months record of the uppermost Hole U1364A thermistor 4 m above seafloor, with an offset in the vertical temperature scale to avoid obscuring the BPR record during that time.
Frontiers in Earth Science | www.frontiersin.org December 2020 | Volume 8 | Article 568566 pressure data and was instrumental in organizing the ONC archiving of the thermistor cable data. JM led the Woods Hole Oceanographic Institution (WHOI) efforts for Keck Foundation funding for the seismic-strain-tilt components of the sensor string and NSF funding for the 2016 installation cruise, and he was chief scientist for that installation cruise. JC was coinvestigator for the WHOI Keck Foundation and NSF funding, participated in the 2016 installation cruise, and was chief scientist on the 2019 cruise to recover the sensor string from Hole U1364A.

FUNDING
KB was supported by NSF grant OCE-1259718 for construction and deployment of the thermistor cable in the hole. Construction of the seismic-strain-tilt instrumentation was supported by a Keck Foundation grant to WHOI, and deployment and recovery of the integrated sensor string was supported by NSF grant OCE-1259243 to JM and JC. Support for the pressure-monitoring instrumentation and 2014 CTD profile was provided by the Geological Survey of Canada and Ocean Networks Canada.

ACKNOWLEDGMENTS
Drilling of Hole U1364A and installation of the ACORK infrastructure was carried out by the Integrated Ocean Drilling Program using D/V JOIDES Resolution. The CTD used in the 2014 depth check was provided by D. Fornari, WHOI. We thank the officers and crews of the JOIDES Resolution, ROPOS, Sikuliaq, and Jason for outstanding support for the at-sea deployments. We thank J. OBrien and K. von der Heydt for engineering support in designing and deploying the WHOI/RSMAS borehole instrument package, and D. Kot, J. Ryder, and I. Kulin for additional technical support during the deployment cruise. An early draft of this manuscript benefitted greatly from comments by H. Villinger.
We thank reviewers C. Ruppel and S. Hickman for careful and constructive comments that prompted significant improvements to this report.