Deep Fault‐Controlled Fluid Flow Driving Shallow Stratigraphically Constrained Gas Hydrate Formation: Urutī Basin, Hikurangi Margin, New Zealand

The Hikurangi Margin east of New Zealand's North Island hosts an extensive gas hydrate province with numerous gas hydrate accumulations related to the faulted structure of the accretionary wedge. One such hydrate feature occurs in a small perched upper‐slope basin known as Urutī Basin. We investigated this hydrate accumulation by combining a long‐offset seismic line (10‐km‐long receiver array) with a grid of high‐resolution seismic lines acquired with a 600‐m‐long hydrophone streamer. The long‐offset data enable quantitative velocity analysis, while the high‐resolution data constrain the three‐dimensional geometry of the hydrate accumulation. The sediments in Urutī Basin dip landward due to ongoing deformation of the accretionary wedge. These strata are clearly imaged in seismic data where they cross a distinct bottom simulating reflection (BSR) that dips counterintuitively in the opposite direction to the regional dip of the seafloor. BSR‐derived heat flow estimates reveal a distinct heat flow anomaly that coincides spatially with the upper extent of a landward‐verging thrust fault. We present a conceptual model of this gas hydrate system that highlights the roles of fault‐controlled fluid flow at depth merging into strata‐controlled fluid flow into the hydrate stability zone. The result is a layer‐constrained accumulation of concentrated gas hydrate in the dipping strata. Our study provides new insight into the interplay between deep faulting, fluid flow and gas hydrate formation within an active accretionary margin.


Gas Hydrate and Seismology
Gas hydrate, a solid clathrate compound of water and gas that is widespread on continental slopes, is a significant component of the global carbon cycle and plays a role in dynamic slope processes (e.g., Sloan, 2003).Gas hydrate stability is primarily governed by (high) pressure and (low) temperature but is also affected by gas composition and pore water salinity.For hydrate to form, gas must be present in concentrations that exceed solubility at the local temperature and pressure.Suitable conditions for gas hydrate occurrence are generally restricted to continental slopes and polar regions (e.g., Kvenvolden et al., 1993).Initial identification of gas hydrate systems is typically undertaken using active source multi-channel reflection seismology.In seismic data, gas hydrate occurrences are generally inferred from observations of a bottom-simulating reflection (BSR) that usually occurs in shallow sediments at the expected depth of the base of gas hydrate stability (BGHS).As such, the BSR is interpreted to represent the transition from free gas in the pore space below the BGHS to gas hydrate within the pore space above it (Shipley et al., 1979).Departures from the expected depth can be used to infer lateral variations in apparent heat flow resulting from fluid advection (Hornbach et al., 2005(Hornbach et al., , 2007)).In this paper, we analyze a uniquely imaged gas hydrate system along a part of the Hikurangi Margin, New Zealand, to better understand the relationship between tectonic structure, stratigraphic controls, and fluid advection that focus on gas hydrate formation.
Many seismic evaluations of gas hydrate accumulations have been undertaken using conventional hydrocarbon exploration techniques, which usually focus on deeper reservoir units using seismic sources with power and frequency properties appropriate for these targets.However, lower-power, higher-frequency sources have also been used for targeted hydrate studies because only shallower penetration depths are necessary (e.g., Haines et al., 2017;Tréhu et al., 2004;Vanneste et al., 2001;Wood et al., 2008).Such investigations have found that gas hydrate systems can have a distinct frequency-dependent response that causes the BSR to decrease in amplitude or disappear altogether with increasing frequency (e.g., Vanneste et al., 2001).This phenomenon is possibly the result of the base of hydrate stability being a transitional rather than a distinct boundary, with the base of hydrate occurring some distance above the actual base of hydrate stability (e.g., Nole et al., 2018) as a result of the processes that are required to convert gas (primarily methane) and pore water into hydrate (Clennell et al., 1999).
The analysis of seismic amplitude variations with offset (AVO) enables better constraints on sedimentary properties (e.g., lithology, petrology, porosity, permeability), pore fluids (e.g., gas, oil, water), and hydrate occurrences within a hydrate setting (e.g., Chen et al., 2007;Ecker et al., 1998).However, analysis of an AVO response requires observations over a sufficiently wide range of source-receiver offsets.Hydrocarbon exploration data are often suitable for such an approach, but short-offset lower-power higher-frequency surveys conducted specifically for shallow-penetration hydrate characterization are often insufficient for such analyses.
Understanding the effects of acquisition parameters such as source configuration and streamer length on seismic resolution is particularly important when analyzing gas hydrate systems because thin stratigraphic (e.g., porous sandstone or tight shale units) or structural (e.g., fault) features affect fluid migration and hydrate accumulation (Boswell et al., 2012;Crutchley et al., 2015;Fraser et al., 2016;Fujii et al., 2015;Navalpakam et al., 2012;Tréhu et al., 2004), thereby significantly impacting the dynamics of gas hydrate systems.Although there have been informative studies in recent years on the effect of survey acquisition parameters on the seismic response from gas hydrate systems (e.g., Crutchley et al., 2023;Haines et al., 2017), detailed analyses of such complementary data sets are not common.

Geological Setting
The Hikurangi Margin is located on an active subduction zone where the Pacific Plate subducts obliquely beneath the Australian Plate off the east coast of New Zealand's North Island (Figure 1).Convergent plate motion decreases from around 45 mm/yr at the northern end of the margin to around 41 mm/yr at the southern end due to counterclockwise rotation of the relative plate motion vectors (Barnes et al., 2010;De Mets et al., 1994).Convergence along the plate interface accommodates roughly 85% of the plate-normal motion (Nicol & Beavan, 2003).The remainder of plate-normal motion and the majority of plate-parallel motion (∼60%) are accommodated within the overriding accretionary wedge through a combination of thrust faulting, strike-slip faulting and vertical clockwise rotation (Nicol & Wallace, 2007).The thrust imbricated frontal wedge of the Hikurangi Margin has formed against a backstop of the Mesozoic Torlesse basement (Barnes et al., 2010;Bland et al., 2015;Lewis & Pettinga, 1993).Thrust faults that extend from the subduction interface into the wedge result in a series of NE-trending anticlinal ridges separated by slope basins that accumulate younger sediment.
The focus of this paper is an intriguing, locally focused gas hydrate system within Urutī Basin (Figures 1 and 2), around 40 km SE of Castlepoint on the coast of the North Island.The distinctive Palliser-Kaiwhata strike-slip fault (Figure 1) cuts obliquely through this region of the wedge, with a local releasing bend on Urutī Ridge (Barnes et al., 1998(Barnes et al., , 2010)).The basement beneath Urutī Basin consists of rocks from the Hikurangi Plateau,  Geochemistry, Geophysics, Geosystems 10.1029/2023GC011144 which is a large Early Cretaceous basaltic oceanic province comprised mostly of volcaniclastic and extrusive volcanic deposits (Mortimer & Parkinson, 1996).Above the basement, the imbricated wedge has a foundation of pre-subduction passive margin Cretaceous-Paleogene rocks in which the décollement is located.These rocks are overlain by Miocene to Recent cover and slope basin sediments.Dredge samples T119 and V479 from the eastern forelimb of Urutī Ridge show that the ridge contains early to middle Pliocene indurated mudstones (Barnes et al., 2010).Urutī Basin sediments were deposited simultaneously with uplift of Urutī Ridge; however, no sediment samples have yet been collected from the basin.

Study Outline
In early 2010, the Urutī Basin gas hydrate system was imaged during the PEG09 petroleum exploration seismic survey (Figures 1 and 2).Urutī Basin contains many seismic characteristics of an active gas hydrate system, such as prominent BSRs and enhanced reflections above and below the BSR that may correspond to hydrate or gas accumulations (Fraser et al., 2016).Fraser et al. (2016) identified a highly reflective possible hydrate-bearing layer within the hydrate stability zone overlying a set of strong reflections beneath the BSR that probably correspond to gas-charged strata.This gas hydrate system was targeted and imaged by the higher-resolution RR1508-HKS02 survey in 2015.
In this paper, we analyze and compare co-located PEG09 and RR1508-HKS02 seismic data over the gas hydrate system to: 1. Interpret the distribution of gas hydrate and underlying gas in this region, 2. Investigate the effects of seismic frequency on the imaging of a focused gas hydrate system to better understand the stratigraphic, structural and dynamic controls that lead to specific configurations of hydrate and gas, and 3. Evaluate enhanced fluid advection through a gas hydrate system using BSR-derived heat flow estimates.
We found that the relationship between stratigraphy, faulting, and fluid flow can provide critical insight into processes controlling gas hydrate formation.

Data Acquisition
The PEG09 survey was conducted in 2009-2010 by Reflect Geophysical with the MV Reflect Resolution under contract from the New Zealand government as part of a regional conventional hydrocarbon exploration study.The acquisition system included three Bolt APG 8500s airguns and a 10,000-m-long 800-channel (12.5 m group interval) hydrophone streamer towed at an approximate depth of 11 m (RPS Energy Pty Ltd, 2010).The shot spacing was 37.5 m for a nominal fold of 133.
In June 2015, the RR1508-HKS02 seismic reflection survey (referred to henceforth as the HKS02 survey) was collected aboard the RV Roger Revelle during voyage RR1508.This survey incorporated a high-resolution shortoffset acquisition system including two 45/90 cu.in.generator-injector (GI) airguns and a 600-m-long 48-channel (12.5 m group interval) hydrophone streamer towed at a depth of approximately 3.5 m with a shot spacing of 25 m for a nominal fold of 12.

Seismic Processing
The seismic data in the vicinity of Urutī Ridge have been processed with the aim of maximizing the resolution of the shallow gas hydrate system while preserving true relative amplitudes.Processing of seismic lines PEG09-25 and HKS02-01 to -06 aimed to minimize the introduction of differences in the stacked data resulting from the greatly different offset ranges in the two data sets.Due to the effects of streamer feathering and limited constraints on the HKS02-01 receiver positions, shot and receiver positions between the two co-located profiles differ slightly (<100 m).To counter this, the CMP positions for the surveys were set to be along straight coincident lines with a common CMP spacing of 6.25 m.Time shifts accounted for the start-of-trace delay and acquisition depth.Based on frequency analysis of the CMP gathers, trapezoidal bandpass filters of 5-10-50-60 and 12-35-160-210 Hz were applied to the PEG09-25 and HKS02 data sets, respectively, to attenuate low-and high-frequency noise.To correct for spherical divergence, a recovery function of t × V RMS (t) 2 was applied to both data sets, where Geochemistry, Geophysics, Geosystems 10.1029/2023GC011144 t is two-way-time and V RMS is the root-mean-square velocity derived from semblance-based velocity analysis of the PEG09-25 data.Spherical divergence gain recovery was preferred over automatic gain corrections as it preserves relative amplitude information.Pre-stack Kirchhoff time migration was then applied to both data sets and a phase shift of 180°was applied to the HKS02 lines to match the PEG09 lines, assuming a zero-degree positive American polarity convention, for example, at the seafloor (Figures 3a,3b,3e,and 3f).

Pre-Stack AVO Inversion
To further characterize the fine-scale structure, we carried out AVO seismic inversion using Hampson-Russell software, which requires a low frequency impedance model (a background model) as input.Typically, background trends of P-and S-wave impedance and density are constrained by well data, but due to the absence of wells in the study area, we derived pseudo-wells at every 100th CMP based on our stacking velocity model.Interval P-wave velocities were calculated using the Dix equation.Densities were determined using Gardner's relationship between P-wave velocity and bulk density (Gardner et al., 1974).
To establish an initial S-wave velocity model, Castagna's "mud rock" equation (Castagna et al., 1985) was used from the seafloor to the BSR, as this region is expected to be dominated by wet clastic units at these shallow depths.However, this empirical relationship underestimates S-wave velocities in gas charged systems, so below the BSR-where the sedimentary units are expected to be more lithified and potentially gas charged-the Krief equation (Krief et al., 1990) was used with coefficients set for a matrix of sandstone with gas as the pore fluid.Note that the empirical rock physics relationships used here have provided acceptable results, but their applicability could be considered in future work (e.g., Avseth et al., 2010).See Supporting Information S1 for further details on the establishment of the initial velocity and density models for the inversion and Fraser (2017) for full details of inversion tests undertaken on these data.
Pre-stack AVO inversion was tested at the pseudo-well locations to validate the input parameters for robustness before inverting the whole seismic section.A typical correlation factor of 0.98 between the real (input data) and synthetic data was achieved for the pseudo wells (Fraser, 2017).P-impedance, Z p , S-impedance, Z s and V P /V S were computed for Urutī Basin data.The V P /V S ratio is more sensitive to a change in fluid type than V P or V S alone.

Heat Flow
To investigate the role of fluid advection through the gas hydrate system, apparent heat flow was derived from the BSR using a simple conductive model (Henrys et al., 2003;Shankar et al., 2010;Townend, 1997).This technique involves calculating the thermal gradient from the observed BSR depth below the seafloor assuming that the BSR marks the free gas to hydrate phase conversion for a particular gas composition.The thermal gradients thus derived can be converted to heat flow by estimating the thermal conductivity.Gridded horizon depths (i.e., seafloor and BSR) were extracted from depth-converted sections to account for the effect of variable velocity (1,500-2,000 m/s) in the shallow seafloor on travel time.
The seafloor temperature, T sf , was derived by extrapolating downward from the lower part of a nearby expendable bathy-thermograph (XBT) drop (drop HKS02-22) collected at the site (Baker, 2016).Temperature at the sea floor as a function of water depth, T sf (d sf ), for these data was thereby estimated to be: where d sf is the seafloor depth.Pressure at the BSR, P, was assumed to be hydrostatic in agreement with other studies in the region (e.g., Henrys et al., 2003;Pecher et al., 2010).
Gas hydrate phase boundary conditions approximated at the BSR were used to determine the temperature at the BSR, T BSR , using the relationship of (Bouriak et al., 2000) for a methane-seawater system.This formulation assumes a gas composition of 100% methane and a standard sea-water salinity of 34-PSU.Gas composition is consistent with previous analyses on the southern Hikurangi Margin (Koch et al., 2016;Schwalenberg, Haeckel, et al., 2010;Schwalenberg, Wood, et al., 2010).
Assuming a constant geothermal gradient between the seafloor and the BSR, the geothermal gradient, G, was calculated using: where d sed is the sediment thickness between the seafloor and the BSR.
To calculate heat flow, the thermal conductivity of the system, κ, was estimated to be 1 W/m°C (cf.Henrys et al., 2003;Townend, 1997) using Hamilton's (1978) porosity-depth relationship for near surface terrigenous sediments and the geometric mean conductivity-porosity relationship of Woodside and Messmer (1961).This allowed heat flow values, H, to be calculated using the simple conductive heat transport relationship: Heat flow values were then gridded using a flex gridding algorithm.Townend (1997) corrected for sedimentation rates and found there to be significant changes in heat flow values, with high sedimentation rates having a dampening effect on uncorrected heat flow values.However, since sedimentation rates are largely unknown in Urutī Basin, heat flow maps were left uncorrected, a stance also taken by Henrys et al. (2003).Additionally, because we are concerned with differences in heat flow within a localized region, a sedimentation rate correction is less important since it should be similar across the region of interest.That said, stratigraphic pinch outs do suggest some differences in sedimentation rates within wider regions undergoing active tectonic changes (see scale of pinch outs in Urutī Basin visible in Figure 3).

Seismic Observations
Seismic lines PEG09-25 and HKS02-01 show similar representations of the gas hydrate system in Urutī Basin (Figures 3a and 3b; position indicated by box in Figure 2).Note several strong reflections, interpreted to correspond to gas charged layers, beneath an inclined BSR that is found beneath a relatively horizontal seafloor.
The strongly reflective unit above the BSR is interpreted to be a layer (or layers) hosting concentrated hydrate.This reflectivity pattern is typical for layer-constrained systems of hydrate overlying free gas observed elsewhere (e.g., Boswell et al., 2012).
Comparison of these two data sets requires a consideration of resolution limitations.In simple terms, the minimum resolvable thickness of a layer is generally accepted to be a quarter of the wavelength (Yilmaz, 1987).This gives a minimum resolvable thickness of 14 m for PEG09 compared to 7 m for HKS02 using a peak frequency of 29.3 Hz for PEG09 and 60.5 Hz for HKS02 and a seismic velocity of 1,680 m/s.Peak frequencies were obtained from frequency spectra extracted from the upper 0.75 s of basin sediment.Lateral resolution is related to the Fresnel Zone, which describes the portion of a reflector where reflected waves interfere constructively (Yilmaz, 1987).For example, for an arbitrary reflection at 1.9 s (the approximate position of the BSR) with an average velocity of 1,530 m/s, the first Fresnel Zone diameter would be 195 m for PEG09 and 136 m for HKS02.The Fresnel Zone increases (and therefore lateral resolution decreases) with increasing depth.Although migration tends to laterally collapse the Fresnel Zone to roughly the dominant wavelength (52 m for PEG09 and 25 m for HKS02) (Stolt & Benson, 1986), lateral resolution may be affected by out-of-plane energy.
The differences in the resolution of the high-and low-frequency data sets enable some specific features to be considered more fully.First, note that in both the low-and high-frequency data, the reflection attributed to a concentrated hydrate unit manifests as a single high-amplitude positive wavelet (Figures 3e and 3f).Since the resolution of the high-frequency data is greater, an interpretation of the hydrate being found in a thinner bed (or set of beds) is therefore possible.Second, consider the effect of resolution on layers within the gas zone.The reflective region below the BSR that is interpreted to be gas charged appears to be more homogeneous in the lowfrequency data than in the high-frequency data where the gas appears to occur in layers that are more distinct and separated from each other.This finding suggests that the interpreted continuity of the BSR depends on the wavelength of a reflecting wave compared to the thickness of layers crossed by the BSR, and is consistent with previous work that has found the observed continuity of a BSR being affected by the frequency content of the seismic data (e.g., Papenberg, 2004;Wood & Gettrust, 2001).Details of the stratigraphy are also much clearer in the high-frequency data; for example, minor erosional surfaces, pinch outs and possible mass flow deposits can be seen better in the high-frequency data (Figure 3b).
Seismic attributes provide further information to assist in the interpretation of the two seismic data sets.Instantaneous amplitude images of the high-and low-frequency data (Figures 3c and 3d) show enhanced amplitudes associated with the horizon inferred to be a porous layer hosting gas hydrate within the hydrate stability zone and charged with free gas below the BSR.As expected, the thickness of high-amplitude layers is generally greater in the low-frequency data, suggesting that bed-thickness tuning effects should be taken into consideration when assessing individual gas zones in the region.Simple wedge models following the method of Hamlyn (2014) suggest that the thickness at which two events become indistinguishable for a single hydrate-bearing layer is around 15.0 m for the PEG09 data set and 4.5 m for the HKS02 data set (Baker, 2016).

Hydrate and Gas Distributions Evaluated by Inversion
AVO inversion of the 2D seismic data from transect PEG09-25 (Fraser, 2017) (Figure 4) provides a means to constrain sedimentary pore space content (i.e., water, gas, hydrate) in a gas hydrate setting (Dutta & Dai, 2009;Fohrmann & Pecher, 2012).However, the lack of good control here requires the introduction of some assumptions concerning the background lithology and stratigraphy based on regional interpretations of sediment types and accumulation patterns (e.g., Barnes et al., 2010;Kroeger et al., 2015).Despite the ambiguities resulting from the lack of definitive lithological controls, the processed seismic images, seismic attributes, and AVO inversion results are consistent with a gas hydrate system associated with a higher porosity and higher permeability sandier unit.See Supporting Information S1 for more details.
Beneath the BSR, an area of low V P /V S ratio suggests the presence of gas (Figure 4d).Within this region is a distinct horizon with relatively high impedance (Figures 4b and 4c), which appears to be stratigraphically continuous with a relatively lower impedance layer above the BSR and is parallel to the overall stratigraphic trend.This horizon may represent a sealing unit that limits fluid transport through it.Above the BSR, high P-and S-wave impedance horizons (Figures 4b and 4c) that correspond to the interpreted hydrate-bearing layer (Figure 4a) may represent stratigraphically controlled hydrate accumulations.These inversion results support the idea of localized fluid flow affecting the depth of the BSR, with at least one low-permeability horizon capped by a layer of lower permeability helping to channel the flow of gas from below, leading to the accumulation of hydrate within the hydrate stability zone.An upward flexure of the BSR, indicative of higher heat flow, is consistent with this interpretation.As the BSR does not appear as a distinct stand-alone reflection in any of the images, it is unlikely to be the result of a wide-spread free gas accumulation, but rather is indicative of the transition from free gas to lack of free gas in permeable strata.

Regional Configuration of the Focused Hydrate Accumulation
The closely spaced network of parallel HKS02 2D seismic lines, spanning a distance of about 7 km along Urutī Ridge and Basin and collected in the area of the hydrate anomaly on PEG09-25 (Figure 5), enable a preliminary assessment of the 3D structure of the feature.The four parallel seismic lines (HKS02-01 to -04) each show: (a) bright reflections below the BSR corresponding to gas charged sediments, (b) an inclined and irregular BSR that in places shallows toward land under a relatively flat seafloor, and (c) a positive polarity reflection (of variable strength) within the hydrate stability zone interpreted to be a hydrate charged layer.
The width and height of the interpreted hydrate accumulation increases from southwest (950 m long and 75 ms high at line HKS02-04) to northeast (2,500 m long and 218 ms high line HKS02-02), suggesting that the northeastern extent of the feature has not been imaged (Figure 6).Also, the maximum uplift of the BSR (measured above the regional BSR to the southeast below the hydrate anomaly) increases from <10 ms on HKS02-04 (Figure 5d) in the southwest to about 100 ms on HKS02-02 in the northeast (Figure 5a).The region of most pronounced BSR uplift imaged in Line HKS02-02 is underlain by a relatively broad (∼2 km wide) region of suppressed amplitudes and acoustic turbidity (Figure 5a).This zone of suppressed amplitudes coincides with the upper extent of the backthrust interpreted by Bland et al. (2015) (Figures 2 and 3).

Heat Flow Implications of the Uplifted BSR
Heat flow estimates based on the regional interpretation of the BSR reveal a distinct heat flow anomaly in the vicinity of Urutī Basin hydrate anomaly (Figure 6).Heat flow is modeled to rise from a local background of about 37 mW/m 2 to a maximum of about 45 mW/m 2 just landward of the intersection of the hydrate anomaly with the regional BSR.The upper extent of the backthrust (Figures 3 and 6) coincides with the position of the heat flow anomaly (Figure 6).Additionally, the most laterally extensive (i.e., broadest) region of the gas hydrate anomaly (imaged on Lines HKS02-01 and HKS02-02) lies updip (i.e., to the southeast) of the most pronounced region of the heat flow anomaly.

Controls on the Depth of the BSR
The PEG09-25 and HKS02 seismic data support the interpretation of upwelling fluids, probably from several km depth in the accretionary wedge, that are channeled upward along deep-seated faults (i.e., the backthrust, Figure 2) to terminate in younger sediments just below the BSR (Figures 2  and 3).The advection of warm fluids at a focused location can result in lateral variability of heat flow, thereby causing a localized uplift in the base of gas hydrate stability (e.g., Zwart et al., 1996).As mentioned previously, the location of the greatest modeled heat flow coincides with the updip termination of the underlying backthrust (Figure 6).
There is a very narrow range of transport rates that can lead to an observable effect on the BSR-too low and there is no visible effect, too high and the BSR is no longer identifiable (i.e., the stability boundary is at or near the seafloor) (Pecher et al., 2010).Using the method of Tréhu et al. (2003), we modeled a range of steady-state flow scenarios that are consistent with our observations at Urutī Basin (Figure 7).In our models, we used the temperature solution of Bredehoeft and Papaopulos (1965), which depends on a given thickness, L, of the upward flow zone (set somewhat arbitrarily to 2,000 m here, corresponding to deep-seated faulting within the accretionary margin) T BGHS is the temperature at the base of the gas hydrate stability zone, T L is the temperature at the base of the upward flow zone, ϕ is the porosity of the sediments, ρ w is the fluid density, c w is the fluid heat capacity, κ m is the thermal conductivity of the matrix, and v z is the vertical flow rate.Our conductive background thermal model uses temperatures at the seafloor and base of the hydrate stability zone of 4 and ∼17°C (based on Bouriak et al. (2000)), respectively.Assuming a thermal gradient of 0.037°C/m, consistent with the regional BSR-derived heat flow, predicts that the base of the hydrate stability zone, z 0 , is ∼350 m below the seafloor, and T L at a depth of 2,350 m below the seafloor will be ∼91°C (Figure 7).Other parameters used for these models were ϕ = 0.33, ρ w = 1,024 kg/m 3 , c w = 4,180 J/kg°C, and κ m = 1 W/m°C, to be consistent with Pecher et al. (2010).
Our mapping at Urutī Basin shows a difference in heat flow of ∼8 mW/m 2 between the peak heat flow and background (Figure 6), and assuming a thermal conductivity of 1 W/m°C (Harris et al., 2019), corresponds to an enhancement of the thermal gradient of ∼0.008°C/m due to the flow.A vertical fluid flow rate of about 1 mm/yr generates a shallow thermal gradient of 0.045°C/m (Figure 7) consistent with the seafloor BSR derived thermal gradient at the fault.Furthermore, assuming a roughly constant temperature of 17°C at the base of hydrate stability that is deflected upward by fluid flow suggests that the BGHS would have risen from 350 to 280 m below the seafloor, which is broadly consistent with the seismic images of the deflected BSR at Urutī Basin (Figure 5).Note that the thickness chosen for the upward flow zone in our model affects the derived flow rates.A shallower source for the fluid would result in cooler temperatures and higher vertical flow rates to generate the same anomaly.
The AVO inversion results show that the shallow stratigraphy in Urutī Basin also contributes to focused fluid flow in this setting.The 3-4 km wide region of low P-wave impedance (Figure 4b), high S-wave impedance (Figure 4c) and low V p /V s (Figure 4d) under the BSR is interpreted to be a zone of gas-charged pore space immediately above the upper termination of the deep-seated fault.However, several landward-dipping strata within this region exhibit P-and S-wave impedances that are markedly higher than the background values.This difference supports geological (i.e., properties related to lithology or porosity) rather than a pore fluid origin for such layers.One probable interpretation is that these high impedance layers correspond to tight shale unitsunder which gas-charged fluids have traveled updip through a porous and permeable layer into the overlying GHSZ.

Constraining a Hydrate Reservoir Unit
Accurate determination of hydrate distributions and concentrations require the identification and characterization of stratigraphic units that can host hydrate or free gas (or other fluids).The quality of such reservoir units depends on their ability to store and deliver the material that they host (Grier & Marschall, 1992).Storage capacity is determined by the volume of interconnected pores (effective porosity) in a reservoir unit, whereas deliverability is a function of a unit's permeability.Within Urutī Basin, several seismically mappable horizons are interpreted to correspond to high-quality reservoir units due to their enhanced amplitudes (and impedance contrasts) within the gas zone beneath the BSR (Figures 4 and 5).By extension, these same units could also be capable of hosting elevated concentrations of gas hydrate within the hydrate stability field above the BSR, with the rationale that significant hydrate growth relies on large interconnected pore spaces (e.g., Torres et al., 2008; Uchida  3).The updip termination of the seismically mapped fault (Figures 2 and 3) is shown by the heavy dashed black line.A high heat flow anomaly closely follows the termination of the fault until the termination reaches a depth of around 2.4 s.The sections of the seismic data shown in Figure 5 are highlighted in white, and the lateral extent of the interpreted concentrated hydrate layer is shown in blue.Map projection coordinates are NZTM2000 in meters.et al., 2004).However, seismic amplitudes alone are not enough to adequately constrain porosity and permeability.Even low gas saturations, if evenly distributed through the pore spaces within a layer, can significantly decrease seismic impedance and enhance reflection amplitudes (e.g., Domenico, 1976).
The distinct high-amplitude feature above the BSR in Urutī Basin is interpreted to result from a high-quality reservoir unit hosting a concentrated gas hydrate accumulation (labeled "hydrate" in Figure 4a).Based on assessments of the stratigraphic setting, the stratigraphic configuration here probably consists of sand-dominated layers interbedded with mud-dominated units within a turbidite sequence (Barnes et al., 2010;Kroeger et al., 2015).This configuration of interbedded sands and shales is interpreted to cause the observed segmentation of the BSR: sandy stratigraphic units promote elevated gas and gas hydrate saturations causing high-amplitude BSR segments, and muddy units inhibit gas and gas hydrate accumulation resulting in weak and absent BSR segments.The high-amplitude feature extends to variable distances above the BSR (Figure 5 and blue highlighted regions in Figure 6).Interpolating between these positions, the reservoir units are observed to cover an area of about 13 km 2 ; however, there are no constraints to the NE or SW, so the full spatial extent of the feature remains poorly constrained.
Estimates of hydrate saturation in Urutī Basin as imaged with seismic line PEG09-25 have been made by Baker (2016) who calculated saturations of 39%-52% assuming clay fractions of 10%-90% (Lee & Deming, 2002).Additional calculations of gas hydrate concentrations in the Pegasus Basin for a range of targets investigated by Fraser et al. (2016) and Crutchley et al. (2016) give values of 49%-56% and 32%-51% for clay fractions of 10%-90%.A more theoretical basin-wide approach to the estimation of hydrate concentrations was made by Kroeger et al. (2015) who predicted that concentrations over 50% would be found in some areas, in agreement with the previous estimates.
High hydrate concentrations such as those calculated above-and manifested as mappable amplitude anomalies within the hydrate stability zone-appear to be restricted to a few layers within the part of Urutī Basin imaged by our data.The seismic data suggest that a combination of factors has led to the establishment of this feature.These include: (a) the presence of a regionally extensive high-quality sand-dominated layers with sufficient porosity and permeability to both facilitate migration of fluids containing natural gas and act as a host for gas hydrate, (b) the superposition of a low-permeability sealing layer that acts as a barrier for vertical fluid flow and facilitates up-dip flow in the permeable strata below, and (c) a source of fluid flow from below providing gas to the system (Figure 8).
In Figure 8, we conceptualize the key processes leading to the concentrated gas hydrate accumulation, starting with fault-controlled fluid flow, and culminating in stratigraphically controlled gas migration into the gas hydrate stability zone.The presence of numerous strata with sufficient porosity is indicated by several highly reflective horizons within the gas charged region.However, the seismic data show that only one main horizon appears to continue as a high-amplitude (flipped polarity) highhydrate-concentration feature within the hydrate stability zone.We interpret that this layer is capped by a sealing unit that acts as a barrier for gas transport from below.The intersection of the good sealing unit with the BGHS results in fluids traveling along the underlying porous and permeable bed.
The high heat flow anomaly based on BSR depth (Figure 6) points to a relationship between an uplifted BSR (Figure 5) and deep-seated faulting below.Furthermore, combining the extent of the high-amplitude hydrate anomaly in seismic data (Figure 5 and indicated by the blue line segments in Figure 6) and the modeled heat flow based on BSR depth shows a link between higher heat flow and the lateral extent of the hydrate anomaly.The greatest heat flow is observed along the northeastern extent of the backthrust in the mapped region, coinciding with the most extensive hydrate anomaly on lines HKS02-01 and HKS02-02.The lowest heat flow is observed along the SW side, corresponding to the shortest hydrate anomaly on line HKS02-04.Since the full regional extent of the upward deflected BSR anomaly has not been mapped by the available seismic data, it is likely that anomalously high heat flow-and therefore hydrate accumulation within the stratigraphically controlled horizon described above-continue to the NE.

Conclusions
An intriguing gas-hydrate anomaly observed in seismic reflection data from Urutī Basin is characterized by a highamplitude reflection extending up into the hydrate stability field from an uplifted BSR.The application of AVO inversion methods supports the presence of concentrated gas hydrate reservoir units within the basin.A grid of high-resolution seismic data shows that the structure of this feature is 3D in nature and significantly affects seafloor heat flow, interpreted to be the result of focused fluid transport from a deep-seated backthrust underlying the basin.Mapping of the extent of the hydrate feature relative to seafloor heat flow (estimated using the depth of the BSR) shows that regions of enhanced heat flow-and therefore enhanced fluid flow-coincide with a more laterally extensive hydrate anomaly.The spatial relationship between the interpreted hydrate anomaly and the underlying gas-charged region suggests that ongoing deep-seated fluid flow plays a significant role in sustaining accumulations of gas hydrate within porous and permeable units.The available data are insufficient to constrain the full spatial extent of the system, so our estimate of the areal extent of the reservoir units is likely to represent a lower bound.

Data Availability Statement
Bathymetric, heat flow, water column temperature, raw seismic data and all other underway data from the RV Roger Revelle voyage RR1508 are available at https://www.marine-geo.org/tools/new_search/search_map.php?&a=1&entry_id=RR1508&output_info_all=on.The calculated heat flow grid shown in Figure 6 and associated ascii files along with all processed seismic data used for this study are available from https://doi.org/10.5281/zenodo.8152964(Gorman, 2023).

Figure 1 .
Figure 1.Urutī Basin study area on Hikurangi Margin southeast of New Zealand's North Island showing seismic data locations over accretionary ridges and basins.Bathymetric data courtesy of Land Information New Zealand cruise TAN1307 (Crutchley et al., 2015).Inset map shows the general plate tectonic and bathymetric setting of the region with continental crust underlying a relatively shallow ocean in light green shades and oceanic crust under a deeper ocean in shades of darker blue.The region of interest highlighted in Figure 6 is indicated by a gray box.Segments of seismic lines shown in Figure 5 are highlighted in yellow.

Figure 2 .
Figure 2. PEG09-25 regional schematic structural interpretation simplified from Bland et al. (2015) highlighting a series of faulted ridges on the Hikurangi accretionary margin.Bottom simulating reflections are indicated by broken yellow lines.The extent of seismic data in Figure 3 within Urutī Basin indicated by pink boxes.Interpreted Urutī Basin backthrust is highlighted in red.

Figure 3 .
Figure 3. Urutī Basin hydrate anomaly in low-(left side) and high-frequency (right side) seismic data.See text for description of specific features labeled in (a).In all images, the top of the deep-seated fault beneath Urutī Basin (Figure 2) is indicated.Horizontal and vertical scales are identical in all sub-figures.(a) PEG09-25 "low frequency" (5-10-50-60 Hz trapezoidal bandpass filter) seismic image.Position of BSR indicated with gray arrows.Position of the hydrate-bearing layer (labeled "hydrate") and gas accumulations are indicated.Note the reversal in polarity of the hydrate-bearing layer as it crosses the BSR: beneath the BSR the reflection is redblack-red (positive-negative-positive), while above the BSR it is black-red-black (negative-positive-negative).(b) HKS02-01 "high frequency" (12-35-160-210 Hz trapezoidal bandpass filter) seismic image.Examples of features imaged more clearly by the high-frequency data are labeled.(c) PEG09-25 instantaneous amplitude.(d) HKS02-01 instantaneous amplitude.(e) PEG09-25 enlargement as indicated in A showing variable area representation of every tenth trace, emphasizing polarity reversal at BSR corresponding to transition from gas fill below to hydrate above.(f) HKS02-01 enlargement as indicated in (b) highlighting same feature as (e).

Figure 4 .
Figure 4. Pre-stack AVO inversion results from PEG09-25 (Boundaries of plots are slightly changed from Figure 3.).(a) Gray-scale image of seismic data from area of interest (cf. Figure 3a) highlighting the hydrate anomaly at Urutī Basin with the background (long-wavelength) seismic velocity (V p ) model derived from NMO analysis overlain in color.Vertical gray lines show the locations of 1D velocity functions that have been interpolated and smoothed to produce the input model for AVO inversion.(b) P-wave impedance (Z p ). (c) S-wave impedance (Z s ).(d) V p /V s .See text for discussion of results.Note that horizontal and vertical scales are identical in all subfigures.

Figure 5 .
Figure 5. 2.5-D characterization of Urutī Basin hydrate anomaly.Seismic images for lines HKS02-02, 01, 03, and 04 (see Figure 1 for locations) show changes in reflective units from line to line across about 7 km of the basin.The broken blue line is the interpreted base of gas hydrate stability (or BSR), mainly identified by the upper termination of high-amplitude (gas-charged) layers.

Figure 6 .
Figure 6.Heat flow calculated over Urutī Basin (colored surface) with structure contours (0.1 s intervals) shown for the underlying fault surface (as interpreted for line PEG09-25 in Figure 3).The updip termination of the seismically mapped fault (Figures 2 and 3) is shown by the heavy dashed black line.A high heat flow anomaly closely follows the termination of the fault until the termination reaches a depth of around 2.4 s.The sections of the seismic data shown in Figure 5 are highlighted in white, and the lateral extent of the interpreted concentrated hydrate layer is shown in blue.Map projection coordinates are NZTM2000 in meters.

Figure 7 .
Figure 7. Modeled temperature profiles for vertical flow rates (v z ) varying from 1 to 4 mm/yr.See the text for details of the model.Increased fluid flow results in the temperature profile bulging upward and a shallowing of the base of gas hydrate stability (as shown by indicative depths z 0 to z 4 in the upper left).

Figure 8 .
Figure 8. Cartoon superimposed on a semi-transparent plot of line HKS02-01 shows a proposed mechanism for the uplifted base of gas hydrate stability (BGHS) and the development of a hydrate saturated sedimentary bed within the hydrate stability field.Fluid flow (blue arrows) is focused along a backthrust underlying Urutī Basin and then directed upward along dipping porous and permeable sedimentary units that cross the BGHS.Interpreted free gas and hydrate pore fill are indicated by light blue dots and pink diamonds, respectively.