Sampling was carried out during July 2019 onboard R/V Tangaroa, as part of the New Zealand HYDEE (Gas Hydrates: Economic Opportunities and Environmental Implications) programme. The study area lies at the southern end of the Hikurangi subduction margin (Figure 1) where the Pacific Plate subducts obliquely beneath the Australian Plate. Measurements were focused on three seeps - Maungaroa, Glendhu Ridge and Uruti Ridge - with data from two additional seep sites–Pahaua Bank and Uruti Ridge West - also considered in the analysis (Turco et al., 2022; see Table 1). Dispersal modelling focused on the main seep site at Maungaroa, where a subsurface layer of gas of >200 m height below the GHSZ has generated doming, fracturing and formation of a conduit that facilitates gas transfer to the seafloor (Crutchley et al., 2021).

A seafloor mooring was deployed in the vicinity of each of the three seeps consecutively to record bottom water physical and biogeochemical properties. The mooring frame was attached to a 460 kg ballast weight and supported a 600 kHz Acoustic Doppler Current Profiler (ADCP) for measuring current direction and velocity. A SBE 37 ODO MicroCAT with conductivity, temperature and optical dissolved O₂ sensors was also deployed on the mooring at ∼15 m above the seafloor. The mooring was deployed for periods of 3–5 days at each site, which either preceded or followed other sampling activities to ensure the mooring and sensors were not compromised. On recovery, the mooring frame was decoupled from the ballast by remote activation of an acoustic release device whilst the ballast remained on the seafloor. The ADCP data was processed and quality controlled to generate a filtered data set with records at 120 s intervals over the lifetime of each respective deployment.

Dissolved methane concentrations were determined by four different approaches: discrete water sample collection in Niskin bottles on a CTD (Conductivity, Temperature and Depth) rosette in either a) standard vertical profiling mode or b) lateral towed mode, with methane concentration also measured by a methane sensor deployed on c) the CTD or d) a Deep Towed Imaging System (DTIS, Bowden et al., 2013). As the dissolved methane plume was largely retained in near-bottom waters, the CTD and DTIS lateral tows focused on the bottom 10 m of water overlying the seafloor, with each sensor tow orientated relative to the seep by reference to a predetermined location and information from preceding tows. Each sensor tow started ∼1 km from the seep and crossed this at the transect mid-point before completion at a similar distance of 1 km. During the CTD sensor tow, Niskin bottle water samples were obtained with increased sampling frequency close to the seep, with the remaining bottles fired in the bottom 200–500 m during retrieval. At least one vertical CTD cast was carried out at a distance of >1 km from the seep at all three sites to characterize background conditions. Bottom waters at Maungaroa received the most intense sampling, with six DTIS tows and 12 CTDs. Stratification was determined using buoyancy frequency squared (N²), which was calculated as: N 2 = g / ρ 0   .   δ ρ / ρ z where g is gravitational acceleration, ρ ₀ is a reference density, and δρ/ρz is the vertical density gradient (Stevens et al., 2012).

Bottom water collected in the Niskin bottles was subsampled into 1-L and 2-L glass bottles on selected CTD tows and profiles at seep and background stations for analysis of dissolved inorganic carbon (DIC) and total alkalinity (TA), respectively. The bottles were overflowed by 100%, with samples preserved by saturated mercuric chloride addition, sealed with parafilm and stored in the dark, with subsequent analysis in the laboratory. DIC concentration was determined coulometrically (Dickson et al., 2007) with an estimated accuracy and precision of ±1 μmol kg^–1 based on repeat analysis of Dickson Certified Reference Materials. TA was determined by closed-cell potentiometric titration (Dickson et al., 2007), using a curve-fitting optimization with least squares analysis of the titration curve, with an estimated accuracy and precision ±2 μmol kg^–1 also based on repeat analysis of Dickson Certified Reference Materials measurements. The DIC and TA data were used to calculate in situ pH on the total scale using in situ temperature, salinity and pressure and Mehrbach dissociation constants (Mehrbach et al., 1973), as refitted by Dickson and Millero (1987).

Stable C isotope ratios of CO₂ were measured using a GasBench II and GC-PAL autosampler coupled to a DELTA V Plus isotope ratio mass spectrometer (IRMS) (all Thermo Fisher Scientific, Bremen, Germany). Daily checks confirmed that IRMS linearity was <0.06‰ volt⁻¹ for a known m/z 44 voltage span, and the standard deviation of δ¹³C values from 10 peaks of working standard CO₂ was always <0.06‰. Values of sample δ¹³C were derived relative to the CO₂ working standard gas introduced directly into the IRMS and corrected for ¹⁷O (Santrock et al., 1985) by Isodat software (version 3.0.94.17; Thermo Fisher Scientific, Bremen, Germany). Isotopic values were reported in δ-notation relative to the Vienna Pee Dee Belemnite scale. This approach was used for both dissolved methane (Section 2.6.1), and also for DIC in the potential oxidation rate measurement (Section 2.7).

Potential methane oxidation rate was measured in selected bottom water samples using ¹³C-CH₄ addition, following the method of Leonte et al. (2017). Seawater was transferred from the Niskin bottle into 6 x 240 mL serum bottles for each depth sample, sealed with butyl rubber septa and aluminium seals, and placed in the dark at ambient bottom water temperature (2–2.5°C). Within three-five hours, four of the six serum bottles were injected with 50 µL ¹³C-CH₄ mixed with 50 µL ambient air from a gas-tight syringe, during which the serum bottle was inverted and the rubber seal perforated by a separate needle to maintain atmospheric pressure in the bottle. The remaining two serum bottles were not amended and provided a time zero for the rate measurement. All bottles were gently agitated for 20 min on a wrist-action shaker, after which the four amended samples were incubated in the dark at ambient bottom water temperature for three-five days. The incubation was terminated at time zero, Day 2 and Day 4 by removal of 2 mL of water; this was injected into a 12 mL gas-tight glass vial sealed with a rubber septum previously flushed with high purity helium and containing 100 µL of phosphoric acid. In the laboratory, acidified samples were equilibrated overnight at 24°C, during which the CO₂ released by acidification of DIC had equilibrated with the headspace. C stable isotope analysis was as described above in Section 2.6. Reference materials IAEA NBS-18 and IAEA-603 were used to normalize all δ¹³C_gas values (Paul et al., 2007). Values of δ¹³C-DIC were calculated from δ¹³C_gas values by accounting for the C stable isotope fractionation associated with the CO₂ gas-aqueous partition (Assayag et al., 2006). Potential methane oxidation rate was determined by the transfer of ¹³C from the methane gas to DIC. Background in situ DIC concentrations were determined using the individual sample CO₂ peak values generated by the GasBench II, calibrated against a subset of samples analysed by coulometry, as detailed in Section 2.5. Rates (nmol l⁻¹ d⁻¹, hereafter identified as nM d⁻¹) were derived using the approach of Leonte et al. (2017), and are regarded as “potential” oxidation, reflecting the large ¹³C-CH₄ addition relative to the background dissolved methane concentration.

Analysis of the integrated discrete water sample and high-resolution sensor data determined that most of the dissolved methane was retained in bottom waters up to 200 m above the seafloor, although this was biased by the large number of sensor measurements obtained within 5 m of the seafloor. Dissolved methane concentration varied over four orders of magnitude, with maximum concentrations exceeding 2 μmol L^-1 (supersaturation of 60,000 relative to atmospheric methane) at the Maungaroa and Uruti Ridge sites (Figure 2, Supplementary Table S1). The lower dissolved methane maximum at Glendhu Ridge may reflect that this site is characterized by a number of smaller dispersed seeps (Turco et al., 2022). All three sites showed methane undersaturation, indicative of methane oxidation, in the mid-water column at depths of 750–1,000 m for the shallower Uruti Ridge site and 1,000–1,500 m for the two deeper sites (Figure 2; Supplementary Table S1). Mid-water column methane concentrations were similar at the background and seep stations at Maungaroa, although elevated methane concentration of ∼20 nmol L^-1 was evident at 1,600–1800 m at the background sites. All three sites exhibited a mean methane concentration of ∼4.5 nmol L⁻¹ between 500 m and the surface, with a maximum of ∼12.5 nmol L⁻¹ at a background site near Maungaroa. The near-surface methane supersaturation of <150% was consistent with previous regional surveys (Law et al., 2010), confirming low regional methane flux to the atmosphere. Although there was no evidence of seep methane reaching the surface, emissions from shallower seeps along the Hikurangi Margin may occur (Higgs et al., 2019). Alternatively, the low methane supersaturation in surface waters may be produced by methane phosphonate or methylotrophic pathways (Karl et al., 2008; Weller et al., 2013).

Regional circulation in the water column overlying the Maungaroa Seep is influenced by eddies and mixing between the northward-flowing Southland Current and the D'Urville Current which flows through the Cook Strait. This generates the Wairarapa Coastal Current (WCC), a relatively cool northward flow along the outer shelf and upper slope (Figure 5D), inshore of the southward-flowing East Cape Current (ECC, Chiswell, 2003). The observed salinity variation between 200 m and 400 m (Figure 5C) likely reflects variability in the relative contributions of the Southland Current, ECC and D’Urville Current between sampling profiles. Below this the Antarctic Intermediate Water (AAIW) is evident as a salinity minimum around 1,000 m depth, with salinity increasing to the seafloor in Pacific Deep Water (PDW).

The mean Buoyancy Frequency squared (N²) profile showed strongest vertical stratification between 250 m and 400 m depth at the base of the WCC (Figure 5B), and then decreased to the core of the AAIW at 1,000 m depth. N² subsequently increased to ∼1,250 m depth, at the interface between the AAIW and PDW, before decreasing again to the seafloor. Comparison with the dissolved methane profile at Maungaroa (Figure 5A; also Figures 4A, B) suggests that water column structure and stability have a limited influence on vertical dissolved methane distribution. The bubble plume extended to a water column depth of ∼1,000 m, shallower than the mid-water column stratification maximum at 1,250 m (compare Figures 5B, E), whereas dissolved methane was largely retained below 1,500 m where there was no corresponding N² maximum (compare Figures 5A, B). The absence of near-seafloor stratification (Figure 5B) is in contrast to Opouawe Bank nearby, where dissolved methane from a shallower seep (depth of ∼1,080 m) was constrained by elevated N² just above the seafloor (Law et al., 2010).

The mooring at Maungaroa was located less than 100 m from the seep centre, with sensor measurements at 15 m above the seafloor. The ADCP recorded relatively low currents (<0.15 m s⁻¹), with north-easterly flow dominating and tidal currents running along the shelf parallel to the coast (Supplementary Figure S1). Fitting of M2 tides confirmed that tidal forcing was the dominant factor, despite an apparent anomaly in flow direction on Day 2. The tidal influence was also apparent in bottom water temperature and salinity, which showed an inverse relationship and variability of <0.2°C and 0.02 ppt, respectively, over the 5-day period. Dissolved O₂ concentration showed minor variation of <2 μmol L^-1 relative to a bottom water mean of 152.5–155 μmol L^-1 (44.8%–45.5% saturation with respect to the atmosphere). This variability primarily reflected temperature-related variation in solubility and there was no temperature or dissolved O₂ signal associated with the seep plume (Supplementary Figure S1), as reported in a previous regional study (Krabbenhoefft et al., 2010).

The modelled methane was scaled by reference to the observed dissolved methane concentration in bottom waters within a 1 km radius around the seep using the Gerris model to estimate the methane source function, which was then converted to a flux by applying this to a source volume. Figures 6A, B show two sensor tows that crossed the Maungaroa seep, with the modelled distribution at the corresponding tidal state shown in Figures 6C, D. Overlaying the transects in Figures 6C, D highlights the overriding influence of the tidal cycle on the spatial distribution of seep-derived methane in bottom waters; despite both transects passing within 50 m of the seep, the modelled methane concentrations and distribution differ markedly due to the tidal state, as shown in the corresponding model simulations.

The footprint of the seep was assessed using different approaches. Previous measurement of water column backscatter indicated a seep footprint of ∼200 m² at Maungaroa (Crutchley et al., 2021), although this is likely an overestimate as the bubble plume will be dispersed over a broader area by bottom currents and advection relative to the actual seafloor release area. An alternative approach to estimating seep spatial extent utilized the decline in meiobenthos abundance from the seep centre to background (D. Leduc, pers. comm.) which provided a horizontal length scale of ∼50 m. As the latter estimate is consistent with the length scale reported for active fluid release at seeps along the Hikurangi Margin (Klaucke et al., 2010), a footprint of 50 m x 50 m was applied with the Gerris simulation, to generate a source function of 4.7 x 10⁻⁶ mol m⁻³ s⁻¹ methane, and so an annual methane emission of 3 x 10⁵ kg at the Maungaroa seep.

Measurements in 15 bottom water samples obtained within 1 km of the three seep sites (Maungaroa, Uruti Ridge and Glendhu Ridge) generated potential methane oxidation rates of 1.79 ± 0.64 nM d⁻¹ (Supplementary Figure S2). Rates were relatively insensitive to bottom water DIC concentration (mean 2058.1 µ mol kg^-1), and also showed no trend with distance from seep location or between the three seeps. In addition, there was no relationship between potential oxidation rate and in situ methane concentration, despite three orders of magnitude variation in the latter (Supplementary Figure S2).

Hydrates release methane with a characteristic low stable isotope (δ¹³C-CH₄) value that results from fractionation during bacterial methanogenesis (Reeburgh, 2007), and consequently elevated methane concentrations in bottom water overlying the Hikurangi Margin seeps were associated with low δ¹³C-CH₄. The three seep sites showed a similar relationship between dissolved methane concentration and δ¹³C-CH₄ in waters between 1,000 m depth and the seafloor (Figure 10A), suggesting that there was no significant differences in either the source or processing of methane in bottom waters at the three sites. However, the more extensive dataset obtained at Maungaroa showed greater variability than the other two sites. Analysis of data from this site for the bottom 200 m of water, where most seep methane was retained (Figures 4, 5), showed a linear relationship between methane concentration and δ¹³C-CH₄ (Figure 10B), as would be expected for lateral mixing of seep-derived methane with background water. However, variability precluded identification of end-member values, and indicated the influence of other factors.

Integration of in situ measurements in the Gerris model and ROMS simulation provided a new approach to estimating methane seep emissions on the Hikurangi Margin. The resulting methane flux for the Maungaroa seep of 3 x 10⁵ kg CH₄ y^-1 is 4%–12% of a previous estimate of 2.5–7.2 × 10⁶ kg CH₄ y^-1 derived from hydroacoustic data (Turco et al., 2022). These different estimates reflect methodological assumptions and uncertainties, as described in Turco et al. (2022) and other studies (Mitchell, et al., 2022). The hydroacoustic calculation likely overestimates the flux, as diffusive loss of methane during bubble rise (McGinnis et al., 2006) was not accounted for, whereas the Gerris model included uncertainty regarding the spatial dimensions of the methane source (see Section 3.3). The dissolved methane measurements in this study may also have included fluid input, although previous regional studies indicate that fluid release is relatively minor in relation to free gas (Luo et al., 2016). Although emissions from seeps are generally consistent in terms of location and magnitude on the Hikurangi Margin (Law et al., 2010), temporal variation in emissions has been reported (Naudts et al., 2010). Consequently, the emission estimates in this study and that of Turco et al. (2022) provide boundaries on annual methane emissions at Maungaroa of 0.3–7.2 x 10⁶ kg CH₄ y^-1. This estimate is high relative to other seeps on Hikurangi Margin, for which annual emissions range from 10² to 10⁴ kg CH₄ y^-1 (Naudts et al., 2010; Luo et al., 2016), but is consistent with a first order estimate derived from hydroacoustic data and residence time for the Opouawe Bank seep located at a shallower depth close to Maungaroa (4.9–9.8 x 10⁵ kg CH₄ y^-1, Law et al., 2010). The revised methane emission for the entire Hikurangi Margin, of 1.1–10.9 x 10⁷ kg CH₄ y^-1, is equivalent to 0.8%–7.9% of New Zealand terrestrial emissions (Ministry for Environment, 2022), and of similar magnitude to the calculated air-sea methane flux for the entire New Zealand Exclusive Economic Zone (4.1 x 10⁷ kg y⁻¹, based upon 130% surface supersaturation; Law et al., 2010). This highlights the critical role of horizontal dilution and the pelagic methane filter in minimizing transfer of seep-derived methane to the atmosphere.

A primary aim of this study was to distinguish the relative contributions of dilution and methanotrophy in determining the distribution and concentration of methane in bottom waters along the Hikurangi Margin. The observed decrease in methane concentration in the vicinity of the seep, and corresponding diffusion and advection in the Gerris model, indicate dilution of two orders of magnitude. The equivalent dilution rate of ∼1 mol m^-3 d⁻¹ is considerably larger than the measured potential oxidation rate of ∼2 x 10⁻⁶ mol m⁻³ d⁻¹ and, consequently, advection and dilution are the dominant determinants of methane concentration for the volume of interest and timescales of days in this study.

However, methanotrophy is an important sink for dissolved methane that not only reduces methane concentrations in the deep ocean but also flux to the atmosphere (Mao et al., 2022). The potential contribution of methanotrophy in this study was examined by two independent techniques: measurement of potential oxidation rate in isotope tracer incubations and interpretation of the dissolved δ¹³C-CH₄ values in bottom waters. The potential oxidation rates were relatively invariant, despite being sampled from within the seep plume through to background locations at all three seep sites (Supplementary Figure S2). This uniformity could reflect tidal retention in bottom waters along the shelf (Figure 9), resulting in a homogenised methanotroph community with uniform oxidation capacity. Although the potential oxidation rate technique relies upon addition of isotopic tracer at concentrations in excess of in situ methane concentrations (Leonte et al., 2018), the observed range and low variability are similar to rates reported at other seeps (5.6 ± 2.3 nM d⁻¹; Leonte et al., 2017; 0–3 nM d⁻¹; Graves et al., 2015; 0–2.5 nM d⁻¹; Weinstein et al., 2016), and marine systems (0–2.5 nM d⁻¹, Steinle et al., 2015; Osudar et al., 2015; 1–10 nM d⁻¹; Mau et al., 2013; 0.01–7.5 nM d⁻¹; Mao et al., 2022). This consistency is surprising, considering the range of environments and ambient methane concentrations, and also the different techniques used to measure methane oxidation. However, it should be noted that rate measurements in this and other studies were predominately carried out at atmospheric pressure, which may alter microbial metabolism and so methane oxidation, relative to that at in situ pressure in the deep ocean (Fortunato et al., 2021).

The absence of a relationship between potential oxidation rate and in situ methane concentration, despite three orders of magnitude variation in the latter (Supplementary Figure S2), has been reported in other studies. However, a first-order relationship has also been reported for some hydrate seeps (Valentine et al., 2001), and hydrothermal plumes (de Angelis et al., 1993), with strong evidence from the Deep Water Horizon wellhead blowout in the Gulf of Mexico after which rates rose by four orders of magnitude (Crespo-Medina et al., 2014). Methanotrophic response is slow (Valentine et al., 2001), with rates increasing ∼1–2 months after methane addition and reaching a maximum after 3–4 months, as confirmed by modelling, oxidation rate and in situ O₂ measurements in the Gulf of Mexico (Kessler et al., 2011; Du and Kessler, 2012; Crespo-Medina et al., 2014). This may reflect slow adjustment of methanotroph metabolism or alternatively adaptation at the community level to elevated methane, which may account for the absence of a relationship between oxidation rate and in situ dissolved methane concentration (Supplementary Figure S2). Alternatively, this lag may reflect other controls of methane oxidation (Steinle et al., 2016), as suggested by shorter response times to methane addition at near-seep sites relative to sites at distance (Chan et al., 2019). As dilution was excluded in the Chan et al. study, the results suggest that near-seep methanotrophs were primed by prior exposure to elevated methane, or alternatively availability of a limiting factor such as trace metals or macronutrients (De Angelis et al., 1993; Beal et al., 2009).

The relatively constant potential oxidation rate across three orders of magnitude in dissolved methane concentration in this study resulted in a range in turnover time of 6–1,000 days, consistent with other systems (7–560 days, De Angelis et al., 1993; 3–1,000 days; Osudar et al., 2015; 500 days to decades; Valentine et al., 2001). However, shorter turnover times have been recorded at seeps where mixing and circulation is constrained (0.3–3.7 days, Leonte et al., 2017; 1.3–24.0 days; Weinstein et al., 2016). Indeed, residence time is an important determinant of methane oxidation rate, with regional geomorphology, seasonal and water mass variation influencing methanotroph distribution, abundance and activity via dilution and current flow (Steinle et al., 2016; Gründger et al., 2021). The Maungaroa seep is located on a relatively flat ridge on the Hikurangi Margin (Turco et al., 2022), without canyons or seafloor morphology that might enhance regional retention and entrainment of bottom water, and so accumulation of dissolved methane, and thus methane oxidation, might be expected to be low. Although this is supported by the linear relationship between methane concentration and δ¹³C-CH₄ in bottom waters, other influences are also apparent (Figure 10B). This may reflect tidal circulation and entrainment of methane in the vicinity of the seep, as visualized in the Gerris and ROMS model simulations (Figures 6 and 9), with fresh seep methane of lower δ¹³C-CH₄ mixing with recirculated “older” methane with higher δ¹³C-CH₄ values due to oxidation. Despite potential methane oxidation rates being low (Supplementary Figure S2) relative to mixing, local retention could enhance methanotroph adaptation to elevated methane, resulting in higher δ¹³C-CH₄ values. Conversely, the variability at Maungaroa in Figure 10B may arise due to methane input from other regional seeps, as illustrated in the ROMS simulation in Figure 9. Methane release from hydrates will not have a uniform regional isotopic value but will vary between sites. In addition, thermogenic methane input has been reported for the Hikurangi Margin (Henrys et al., 2009), which could be a source of methane with higher δ¹³C-CH₄. Although tidal circulation of bottom waters may increase methane oxidation and mixing with methane from other sources this also confounds determination of their relative contribution to local methane dynamics.

Methane oxidation provides a source of energy and C for methanotroph growth (Hanson and Hanson, 1996), and ultimately produces CO₂ whilst consuming O₂. This is reflected in the elevated O₂ consumption in seafloor seep communities relative to non-seep benthos (Boetius and Wenzhöfer, 2013), that arises directly from methanotrophy and indirectly from chemosynthetically-supported benthic communities. However, there was no evidence of the Hikurangi Margin seeps influencing bottom water O₂ concentrations in either the mooring data (Supplementary Figure S1) or CTD profiles (data not shown), and the measured potential methane oxidation rates would have negligible effect on bottom water O₂. If the total methane release by the five seeps was completely oxidized within the modeled plume volume (90% of dye tracer, ∼1,000 km³) within 1 year, the decrease in dissolved O₂ concentration would be <0.3%. This assumes utilization of 2 mol of O₂ during methanotrophy and so conversion of all methane-fueled biomass to CO₂ (Chan et al., 2019), with resulting CO₂ production of <0.2 μmol L^-1 and a corresponding decrease in pH of <0.002 relative to a background of 7.87. This minor influence is consistent with observations at the Hudson Canyon seeps, where aerobic methane oxidation accounted for only 0.3% of observed changes in DIC (Garcia-Tigreros and Kessler, 2018). A previous estimate based upon all Hikurangi Margin seeps suggested a greater impact of methane oxidation on dissolved O₂ and the carbonate system (Turco et al., 2022), but this was derived using a higher seep methane flux and also did not consider dilution and dispersion. Although the current assessment indicates negligible impact of seep methane on bottom water biogeochemistry it does not account for superimposition of plumes from different seeps along the Hikurangi Margin, as shown in Figure 9, which may result in localized impacts. For example, oxidation of seep methane in overlapping plumes could generate localized decreases of 0.53% in dissolved O₂ concentration and 0.005 in pH. These localized methanotrophy-driven decreases would compound the current low dissolved O₂ saturation in deep water along the Hikurangi Margin, and contribute to ongoing deoxygenation and carbonate saturation horizon shoaling (Feely et al., 2012; Oschlies et al., 2018).

This study highlights the value of assessing the influence of methane seeps on bottom waters by integrating observations within a hydrodynamic framework. Although bottom water currents confound characterization of dissolved methane distribution around a seep, as illustrated in Figure 6, tidal dispersion and dilution was accounted for in a Gerris model. The ROMS model subsequently simulated methane concentration and distribution using the source function and location generated by the Gerris model but did not include loss to methane oxidation. Although rate measurements indicate that this process is insignificant relative to dispersion, further model development should incorporate methane oxidation, informed by measurements along the Hikurangi Margin to investigate the apparent uniformity in rate and also establish the relative contribution of oxidation to bottom water methane dynamics. The superimposition of multiple seeps could be further examined by applying the ROMS framework to the full suite of seeps along the Hikurangi Margin (Watson et al., 2020). As with many global studies, methane release from seeps is not apparent in surface waters along the Hikurangi Margin, except for some shallow seeps (Higgs et al., 2019); however, there is the potential for transfer via vertical mixing associated with geomorphological features and eddies. Application of the ROMS model framework to the broader shelf region could determine whether this represents a “shortcut” for seep methane transfer to surface waters and flux to the atmosphere along the Hikurangi Margin. The ROMS simulation also suggests that tidal retention and superimposition of plumes may facilitate connectivity between seeps, which should be examined with respect to larval dispersion and connectivity between seep populations (Levin et al., 2016) along the broader Hikurangi Margin. In addition, integration of regional hydrodynamic models with Earth System Models (Rickard et al., 2023) could examine future climate-driven changes in bottom water and the associated biogeochemical impacts along the Hikurangi Margin.