Estimates of Methane Release From Gas Seeps at the Southern Hikurangi Margin, New Zealand

Introduction

Methane forms in marine sediments that are rich in organic matter through either microbial methanogenesis or thermogenic processes (Schoell, 1988). Methane formation occurs at different depths in the subsurface but, because of the buoyancy of the gas, it migrates upwards through pathways that include permeable carrier sedimentary units, faults, or densely fractured regions (Cook and Malinverno, 2013; Crutchley et al., 2015; Nole et al., 2016; Hillman et al., 2017; Hoffmann et al., 2019; Hillman et al., 2020). When these pathways connect to the surface, gas bubbles escape the seafloor as gas seeps, which can range from diffusive sporadic and localized emanations of bubbles to widespread, vigorous gas seeps, occurring in different geological contexts, from the coastal environments to deep ocean regions (Judd, 2004; Duarte et al., 2007; Watson et al., 2020). In shallow waters—up to 800 mbsl, the gas bubbles can reach the sea surface and release greenhouse gases into the atmosphere (Schmale et al., 2005), having direct implications for climate. In deeper seas, however, most of the gas that is released at the seabed dissolves into the ocean without reaching the sea surface, causing localized methane-induced seawater acidification over long time scales (Law et al., 2010; Biastoch et al., 2011; Garcia-Tigreros et al., 2021).

Seismic and acoustic methods are useful tools to identify, map and characterize free gas accumulations in the subsurface (e.g., Judd and Hovland, 1992; Kim et al., 2020), evidence for past and present seepage at the seafloor (e.g., Stott et al., 2019) and gas release into the water column (e.g., Colbo et al., 2014; Böttner et al., 2020). In the subsurface, the presence of gas in the pore-space significantly affects the elastic properties of the bulk sediment, primarily by reducing the bulk seismic velocities and generating a contrast in acoustic impedance. The detection, mapping and characterization of subsurface gas reservoirs are parts of the traditional hydrocarbon exploration workflow (Yilmaz, 2001). Low seismic velocities, often associated with fluid migration through geological structures such as chimneys, conduits, and faults, can be an indicator of free gas in the sediments.

In the ocean, columns of rising gas bubbles (also called flares for the typical shape they assume in hydroacoustic imaging) are less dense than seawater and thus represent strong acoustic reflectors, as they generate sharp localized changes in the acoustic impedance of the water column. Calculating methane fluxes at a seep site is challenging, as it requires ground-truth information about bubble size distribution, chemical composition, density, bubble coating, and ascending speed (Leblond et al., 2014). Typically, gas bubbles released at the seafloor tend to have a radius in the range of 1–15 mm (Veloso et al., 2015), resonating at frequencies from a few hundred Hz to ∼12 kHz, depending on depth and size (Weidner et al., 2019). Single-frequency sonar systems have been successfully used for the identification of seep sites and for water column imaging for decades (e.g., Merewether et al., 1985; Hornafius et al., 1999; Nikolovska et al., 2008). However, these instruments cannot be used to determine the size distribution of the bubble population.

A common approach for the estimation of gas fluxes is coupling acoustic imaging of the gas bubbles with optical point-source measurements from towed camera systems (Higgs et al., 2019), remotely operated vehicles (Naudts et al., 2010), bubble observation modules (Bayrakci et al., 2014) or bubble traps (Römer et al., 2012). Although these point-source measurements provide the most accurate observation of bubble parameters, they require long deployment durations and a restricted field of view of less than ∼15 m. Moreover, they are also limited to measurements at the seafloor, and cannot provide a way to track the changes in bubble size distribution as they rise through the water column. Broadband hydroacoustic methods provide a more efficient tool to directly estimate bubble parameters by insonifying large areas of the oceans using a range of frequencies (e.g., Veloso et al., 2015; Colbo et al., 2014; Dupré et al., 2015; von Deimling et al., 2011; Li et al., 2020).

The highest concentration of cold seep sites worldwide has been observed along convergent margins, (Saffer and Tobin, 2011; Suess, 2020; Watson et al., 2020). In these regions, gas flares observed at the seabed and other shallow gas migration features are often connected to subsurface methane reservoirs through gas conduits (e.g., Meldahl et al., 2001; Petersen et al., 2010; Krabbenhoeft et al., 2013; Crutchley et al., 2021). On the seafloor, gas seeps are the most common manifestations of ongoing subsurface fluid flow (Judd and Hovland, 2009). The gases that are expelled from gas seeps on continental margins are primarily composed of methane, leaving major questions open on: 1) the amount of methane reaching the ocean surface (McGinnis et al., 2006; Shakhova et al., 2010; Fu et al., 2020), 2) the connectivity of seeps to deeper hydrocarbon systems (Crutchley et al., 2021), 3) the role of gas hydrate dissociation (Reagan et al., 2011), 4) how gas flux rates change over time and the potential influence of seismicity on subsurface fluid flow (Bassett et al., 2014; Bonini, 2019; Legrand et al., 2021). The southern Hikurangi Margin, off the North Island of Aotearoa/New Zealand, reveals evidence of widespread methane seepage (Greinert et al., 2001; Barnes et al., 2010; Watson et al., 2020). The accretionary wedge here consists of a series of thrust-related ridges striking NE-SW, composed of compressed and deformed sediments, probably turbidites and ancient trench-fill deposits (Kroeger et al., 2015; Lewis et al., 1998). Seismic data show evidence of concentrated gas hydrate accumulations in the core regions of many thrust ridges in this region, indicating the preferential migration of free gas along permeable strata towards the core of anticlinal structures (Crutchley et al., 2019; Turco et al., 2020; Wang et al., 2017a; Schwalenberg et al., 2010; Kroeger et al., 2021). Estimation of gas fluxes at gas seeps on the Hikurangi Margin has so far relied on in situ optical measurements of gas bubbles (Naudts et al., 2010; Higgs et al., 2019) combined with single-beam acoustic data.

In this paper, we use a combination of seismic and acoustic data to characterize five known sites of active gas seepage on the southern Hikurangi Margin. The interpretation of seismic amplitudes provides a means to identify regions of free gas accumulations within the sediments. Qualitative analysis of multibeam data collected over 3 years allows the spatial extent of the region of active venting at the seafloor to be mapped, while also imaging the acoustic flares in the water column and analyzing their variability over time. Split-beam echo sounder data are used to extract important parameters for the quantification of gas bubbles in the water from the backscattered acoustic energy, which are then used to calculate gas fluxes at the cold seep locations.

Data and Methods

Acoustic Data Processing

The identification and mapping of gas seeps in deep waters was achieved through the analysis of acoustic data (Figure 1). Bathymetric and water column data were acquired during three scientific voyages onboard the R/V Tangaroa (Figure 2): TAN1808 (September-October 2018), TAN1904 (July 2019) and TAN2012 (November 2020).

FIGURE 1. Illustration of the workflow to estimate gas fluxes from MBES (EM302) and single-beam (EK60) backscatter data. Approach (1) is based on the manual curve fitting of the normalized frequency response of the 18 and 38 kHz channels of the EK60 data, while approach (2) is based on a linear inversion of the non-normalized frequency response. S_v, mean volume backscattering strength computed for each cell; $\overline{\mathit{TS}}$, mean target strength computed for the whole acoustic flare; ρ, density of bubbles; BSD, bubble size distribution; v, bubble rising velocity; A, spatial extent of acoustic anomaly.

FIGURE 2. (A) Overview of the Hikurangi subduction Margin of New Zealand (dashed black line) and the study area (yellow rectangle). (B) expansion of the study area on the southern Hikurangi Margin: white, cyan and magenta dots are gas flares identified from the MBES data from the three R/V Tangaroa voyages (TAN1808, TAN1904, TAN2012). The five target areas are mentioned in the text. Map coordinates are in metres of New Zealand Transverse Mercator 2000 (New Zealand Geodetic Datum 2000).

Swath bathymetry and acoustic backscatter of mid-water reflectors were collected with a hull-mounted Kongsberg EM302 multibeam echo-sounder during the three voyages. The EM302 echo sounder operates at a nominal frequency of 30 kHz and with a swath of ∼120°. The use of the multibeam data was twofold: 1) to accurately locate gas seeps on the seafloor and 2) to calculate the area of seepage at the seafloor. For the former objective, the data were processed using the National Institute of Water and Atmospheric Research (NIWA) custom-built software Espresso with the following steps: seafloor detection filtering, removal of the outermost noisy beams (>45°), removal of bad pings, filtering side lobe artefacts and muting the first 5 m of data above the automatically picked seafloor, to avoid misinterpreting the smearing of the beams at the seafloor as gas bubbles (Schimel et al., 2020). The correct pinpoint (as well as could be determined) of seepage at the seafloor was facilitated by the fan visualization of MBES data (Figure 3). To calculate the total area of gas seepage in proximity of the seafloor, the processed data were vertically summed over a window between 5 and 20 m above the seabed - a process known as echo integration (MacLennan et al., 2002). The output of this process is a georeferenced image of volume backscatter intensity with a horizontal spatial resolution of 15 m × 15 m, that allowed mapping the spatial extent of the acoustic anomaly in the water column.

FIGURE 3. Backscatter intensity images of the Glendhu seep field (see Figure 2 for location) from multibeam data. (A) echo integrated map (resolution 20 × 20 m); (B) range stacked view along the profile shown in a); (C) fan view at the location indicated by the dashed black line in (A).

A suite of 5 Simrad EK60 echo sounders were used to obtain calibrated acoustic measurements of the water column during TAN1904 and TAN2012 voyages. The data were acquired over the five targeted areas based on existing multibeam coverage. These split-beam systems were calibrated using a standard 38.1 mm tungsten sphere hung under the vessel, following standard procedures (Demer et al., 2015). Given the relatively great water depth of most targeted flares (>1,000 mbsl), only the 18 and the 38 kHz echo sounders had sufficient range to image the full vertical extent of the gas bubbles. The two-way beam angle of the EK60 echo sounder is 11° for the 18 kHz system and 7° for the 38 kHz one. We extracted and processed the data recorded for the targeted gas flares using the open-source software package ESP3 (Ladroit et al., 2020). The data were processed to only extract the acoustic signal associated with gas venting. The processing included seafloor echo detection and removal, bad ping removal and de-noising (De Robertis and Higginbottom, 2007). Once acoustic flares were identified and extracted, we carried out frequency analysis on the pre-processed split-beam data and compared the frequency response to theoretical bubble backscatter models to estimate the bubble size distribution (BSD) of the entire flare (Figure 1). Finally, we echo integrated the processed 18 kHz data using cells 25 m high and 10 m wide, in order to retrieve a mean volume backscattering strength ($\overline{S_v}$) response for each cell.

Estimation of Gas Fluxes

Existing theoretical models to predict the acoustic backscattering cross-section (${\sigma }_{bs})$ of underwater bubbles differ in terms of parameters such the bubble size distribution (BSD), bubble morphology, the frequency range of the acoustic signal and the environment characteristics. For the scope of this work, we adopt the model proposed by Ainslie and Leighton (2009). For a more complete overview of the definitions and terminology used in underwater acoustics and adopted in this work the reader is referred to, e.g., MacLennan et al. (2002).

The workflow to estimate gas fluxes at the seep locations is illustrated in Figure 1. To estimate the BSD and density, we followed two approaches, one based on the normalized frequency response of the 18 and 38 kHz channels of the EK60 data, and one based on the non-normalized frequency response, similarly to the approach used by Li et al. (2020).

In the first approach (1), gas flares were isolated using the 18 kHz data, and they were echo integrated using a variable cell size ranging from 25–50 m in height and 5–10 m in width. The echo integration process yields a mean volume backscattering strength ($\overline{S_v}$) for each cell of the acoustic flare. The frequency analysis of the gas plume from both the 18 and 38 kHz channels, normalized relatively to the 38 kHz $S_v$ values, provides frequency-dependent $S_v$ values that the theoretical backscatter model curves aim to match. The curve-matching procedure was achieved through manual tuning of the mean bubble radius and bubble size distribution. When the match is considered satisfactory, i.e., when the residuals between the observed and the theoretical $S_v$ are below a chosen threshold for each water depth, the resulting BSD is considered representative of the entire flare, and the mean backscattering cross-section is computed following:

where $BSD$ is the bubble size distribution, and $a$ the bubble radius. ${\sigma }_{bs}$ is given byAinslie and Leighton (2009):

where $\omega$ is the frequency of the active source of the acoustic signal, ${\omega }_1$ is the resonance frequency of the bubble and $\delta$ is a dimensionless frequency-dependant damping constant. To retrieve the density of bubbles per volume unit, the mean volume backscattering strength $\overline{S_v}$ value of each cell of the echo integrated flare is compared to the representative mean target strength ($\overline{TS}$) value of the whole flare. Eq. 3 relates the mean backscattering cross-section representative of the entire flare $\overline{{\sigma }_{bs}}$ (in m⁻¹, or expressed in decibels by $\overline{TS}$) to the mean volume backscattering coefficient of a given cell of the echogram $\overline{s_v}$ (m², or $\overline{S_v}$ in decibels), following:

Then, $\rho$ is averaged across cells located at the same water depth, hence a mean density of bubbles is estimated for each horizontal slice of the gas flare ($\overline{\rho }$).

The second approach (2) to retrieve BSD and $\overline{\rho }$ for each flare consists of a linear inversion technique that aims to minimize the difference between the modelled and the observed non-normalized volume scattering strengths at each depth of the acoustic flare. The forward modelled data are computed following the theoretical backscatter model by Ainslie and Leighton (2009), and the optimization is achieved through a conjugate gradient technique. The model parameters that are inverted are bubble size, bubble size distribution parameters (for a fixed distribution type) and density of bubbles per volume unit. These parameters are initialized using the results from the first approach and updated iteratively during the inversion. 100 iterations were performed on each flare, and the inversion was carried out twice, once assuming a log-normal type bubble size distribution (Veloso et al., 2015), and once assuming a Weibull type distribution (Dey and Kundu, 2012).

Once the BSD and $\overline{\rho }$ are known, the gas flux is estimated for the entire flare (which could represent a field of seeps) following Eq. 4:

where $A$ is the spatial extent of the acoustic anomaly in the proximity of the seafloor as interpreted from the processed MBES data, $\overline{\rho }$ is the mean density of bubbles in the water in m⁻³, $\text{v}\left(\text{a}\right)$ is the mean rising velocity of the gas bubbles, $BSD$ is the bubble size distribution and $a$ is the bubble radius. The bubble rise velocity valuesare based on the work by Leifer and Patro (2002) and consider the two endmembers clean and coated bubbles, where the latter model represents gas bubbles which are coated with oil or hydrate (Fu et al., 2020).

In the next section, the estimated gas fluxes are presented as ranges of values. The major source of variability in the flux estimations comes from the use of coated versus clean bubbles models: because clean bubbles rise faster than coated bubbles, changes in $\text{v}$ have a linear effect on the fluxes (Eq. 4). Hydrate-coated bubbles have been observed at some of the seep sites analyzed in this paper (Sarah Seabrook, personal communication). However, due to the scarcity of observational data, we use both models to provide a range of possible gas flux estimates.The secondary source of uncertainty in the estimations is related to the approach used for the calculation of the $BSD$ and density (Figure 1). These parameters impact both the mean rising velocity–which depends on the bubble size, and the volume of gas calculated at the seafloor.

Because of the lack of in situ chemical measurements at the locations of seepage, we assume that 100% of the gas released at the seeps is CH₄.

Seismic Data

High resolution seismic reflection data were acquired during the TAN1808 research cruise (Figure 2). A GI gun and a 600 m long streamer of 48 channels were used for the acquisition (Crutchley et al., 2018). The seismic data processing is described in detail by Turco et al. (2020) and included geometry application, Butterworth filtering (corner frequencies of 7, 14, 150, and 200 Hz), constant high-dip noise removal through FK filtering, corrections for spherical divergence correction (1,500 m/s), CDP sorting and NMO correction (1,500 m/s), stacking and post-stack Kirchhoff time migration. Due to the presence of gas hydrates within the sediments, we used a seismic velocity of 1700 m/s, based on the velocity analysis by Turco et al. (2020) to depth convert the seismic sections. The dominant frequency of the processed seismic data is 95 Hz, and the vertical resolution is approximately 4–5 m.

Results

Identification of Gas Seeps

Gas seep sites are identified in the multibeam data by anomalously high acoustic backscatter in the water column with respect to the surrounding region. High backscatter values in the water could also indicate the presence of schools of fish, thermo-cline layering or artefacts. Given the ambiguity in interpreting vertically summed (echo integrated) backscatter intensity maps, we analysed horizontally stacked sections and fan-view images of backscatter intensity (Figures 3B,C, respectively) in the vicinity of the acoustic anomalies, to confidently interpret gas flares where regions of high backscatter intensity propagate from the seafloor upwards, as expected from a rising aggregate of gas bubbles. We analysed three datasets from different voyages that surveyed the same target areas (Figure 2). This approach ensures that we accurately pinpoint the location of gas venting at the seafloor and gives a temporal dimension to the study. In the study area, we identified a total of 129 individual gas flares: 46 from TAN1808, 32 from TAN1904 and 53 from TAN2012 datasets (Table 1). Most of these flares are located approximately at the same point on the seafloor in the three datasets; however, the difference in data quality and acquisition parameters does not allow a more detailed comparison between the three datasets. It is important to note that the flares observed in the acoustic data are presumably formed by multiple bubble outlets sited in an area that is smaller than the insonified seabed area. The lateral resolution of the MBES data at the seafloor depends on several factors such as beamwidth, water depth, survey speed and swath coverage. For our study area, the lateral resolution varies between 25 and 50 m. The five regions of focused gas seepage are: Urutī Ridge, Urutī Ridge West, Pahaua Bank, Glendhu Ridge and Mungaroa Ridge (Figure 4). Mungaroa Ridge is an informal name, which has not been officially gazetted by the New Zealand Graphic Board—Ngā Pou Taunaha o Aotearoa. The shallowest seeps occur at Urutī Ridge West at ∼1,100 mbsl, while the deepest one is the Honeycomb Ridge seep, located in ∼2,400 m water depth. There is no acoustic evidence that any of the analysed gas flares reaches the sea surface in the study area.

TABLE 1. Overview of seeps observations based on MBES data.

FIGURE 4. Distribution of the gas flares identified in the five target areas: (A) Mungaroa Ridge; (B) Glendhu Ridge; (C) Urutī Ridge; (D) Urutī Ridge West; (E) Pahaua Bank. The approximate direction of plate convergence is extracted from the MORVEL online tool (Argus et al., 2011). MCS, multi-channel seismic data.

Gas Fluxes and Seismic Observations

In this section, we present the results of gas flux estimations for the five target areas (Figure 5) and analyse the local geological structure of these sites.

FIGURE 5. Results of methane flux estimations for the five seep locations. (A–E) show the variations of the mean bubble radius versus depth; (F–J) show the mean density of bubbles and (K–O) show the mean methane fluxes calculated at depth with the clean and coated ascending velocity models. 2 W: linear inversion with Weibull distribution; 2 L: linear inversion with log-normal distribution; 1 W: manual curve fitting with Weibull distribution.

At each seep site, we selected one representative flare from the SBES data (shown in Figures 7–10). The selection was driven by data quality, representativeness of the flare for the entire seep field, and vicinity to the location of the seismic line. The selected flares were used to estimate bubble size parameters and bubble density, which were considered to be representative of all the seeps located in the same field.

The flux estimates provided in the following sections represent an average of the linear inversion method (Approach 2), while the comparison of the results of the manual curve fitting (Approach 1) is shown in Figure 5 and in Table 2, together with the details of the parametrization of gas flares for each of the target sites.

TABLE 2. Details of bubble size parameters, density of bubbles, mean gas fluxes and seep water depth for the five target areas.

The average thickness of the gas hydrate stability zone (GHSZ) varies according to the water depth and the geological context, from ∼360 m at Urutī Ridge (1,200 mbsl) to ∼630 m at Mungaroa Ridge (2,100 mbsl). Despite being visible throughout the five study areas, bottom-simulating reflections (BSR) associated with gas hydrate occurrence are discontinuous and cannot be observed directly below the locations of gas seepage. In Figures 7–11, the vertical and horizontal scales of the EK60 echograms showing the gas flares are approximately equal to the scales of the seismic image. It is important to note that the acoustic images of the water column show the apparent resolution of the EK60 data: the true horizontal resolution depends on the beamwidth and varies with depth following:

Where r is the horizontal resolution, r^∗ is the apparent resolution, d is the water depth and θ is the beam aperture angle (Figure 6). Seismic velocity of 1700 m/s (Turco et al., 2020) was used to depth convert the seismic sections, whereas the echograms were depth-converted on-the-fly during data acquisition using sound velocity profiles. The dominant frequency of the processed seismic data is 95 Hz, and the vertical resolution is approximately 4–5 m. The acquisition and processing parameters of the seismic data are provided by Crutchley et al. (2018) and Turco et al. (2020), respectively.

FIGURE 6. (A) Graphical illustration of the lateral resolution of the acoustic images of the single-beam (EK60) data. (B) Plot showing the dependence of the lateral resolution (L) on water depth and beam angle.

Urutī Ridge

The main Urutī Ridge seep field is located slightly seaward of the bathymetric high of the anticlinal ridge. The seepage occurs over ∼4 km² of the seafloor, in water depths from 1,175 m to 1,300 m, and tens of flares can be identified from the acoustic data (17 flares identified in the 2018 datasets, 14 flares in 2019 and 9 in 2020). The same flares were imaged in all surveys. Most acoustic flares reach ∼700 m water depth, and the total area of high acoustic anomaly at 20 m above the seafloor measures 0.43 km². The easternmost flares imaged in the TAN2012 dataset seem to be aligned roughly NS, which is sub-perpendicular to the direction of plate convergence in this part of the margin. The main flare was selected from the total population of flares of this region to estimate the density of bubbles per cubic metre. Gas flux estimates for the entire seep field at Urutī Ridge range between 1.28 and 3.55 × 10⁶ kg/year. The seismic profile shown in Figure 7 runs perpendicular to the strike of the main anticlinal structure, and it crosses the seabed location of two major gas flares used for flux estimations. A broad extent of ∼1 km of the shallow subsurface shows high-amplitude negative polarity reflections that reveal the presence of free gas in the sediments. The sedimentary sequence below this region is characterized by a general decrease in seismic amplitudes (seismic blanking) and disrupted reflections. The blanking zone in the overlying stratigraphic sequence is bounded in depth by a seismic unconformity that marks the top of a highly reflective unit of steeply dipping strata that form the seaward limb of the Urutī Ridge anticline. The amplitude of the BSR is higher to the NW and to the SE of the flare site, it fades out in the core region of the anticline, and it is not observed in the region of enhanced reflectivity on the seaward limb of the anticline.

FIGURE 7. Overview of the Urutī Ridge seep site: the interpreted TAN1808-22 seismic profile is shown in the main panel. The bold black arrows represent the direction of fluid flow as interpreted from the seismic data and explained in the text. BSR, bottom simulating reflection. (A) Map view of the Urutī Ridge seeps (red box in Figure 4C), showing the acoustic backscatter anomaly in the echo integrated MBES (EM302) data in proximity of the seafloor. The location of the TAN1808-22 seismic line and of the single-beam data (EK60) are indicated by the black and the red lines, respectively. The black dots represent the location of the main gas flares. The hydroacoustic data were collected during the TAN1904 voyage. (B) Echogram of two gas flares as imaged by the 18 kHz channel in the single-beam data. The horizontal scale is the same as the seismic panel. (C) Expanded view of the seismic data showing the shallow region beneath the cold seeps.

Glendhu Ridge

Glendhu Ridge is a thrust-related elongated structural feature with four-way closure that lies close to the present-day deformation front. The anticlinal structure of the ridge is imaged in the seismic profiles and has been analysed in detail by Turco et al. (2020). There is no BSR below the seep location at the top of the ridge, similar to what is observed at Urutī Ridge (Figure 7). The main gas venting field is located right on the bathymetric crest of the ridge, at a water depth of about 2000 m, where 6–8 main acoustic flares can be identified from the multibeam data. The seeps are roughly aligned ENE-WSW, parallel to the long-axis of the four-way closure and sub-parallel to the vector of plate convergence (Figure 4B). For the parametrization of this gas seep field, we used the main acoustic flare visible in Figure 8, which rises from the seafloor for roughly 1,200 m reaching a depth of ∼780 mbsl. The area of acoustic anomaly at this site is 0.17 km², yielding total gas flux estimates of 1.57 and 4.34 × 10⁶ kg/year, considering coated and clean bubbles, respectively.

FIGURE 8. Overview of the Pahaua Bank seep site: the interpreted TAN1808-44 seismic profile is shown in the main panel. The bold black arrows represent the direction of fluid flow as interpreted from the seismic data and explained in the text. BSR: bottom simulating reflection. (A) Map view of the Pahaua Bank seeps (red box in Figure 4E), showing the acoustic backscatter anomaly in the echo integrated MBES (EM302) data in proximity of the seafloor. The location of the TAN1808-44 seismic line and of the single-beam data (EK60) are indicated by the black and the red lines, respectively. The black dots show the locations of the main gas flares. The hydroacoustic data were collected during the TAN1904 voyage. (B) Echogram of a gas flare as imaged by the 18 kHz channel in the single-beam data. The horizontal scale is the same as the seismic panel. (C) Expanded view of the seismic data showing the shallow region beneath the cold seeps, where free gas accumulation is inferred by the negative polarity reflection.

Pahaua Bank

Pahaua Bank is a submarine ridge located on the mid-slope portion of the accretionary wedge, at water depths of 1,450–1,570 m. There are two regions of gas seepage at the seafloor: the northernmost group of gas seeps includes at least seven distinct flares aligned NNW-SSE, perpendicular to the direction of plate convergence (Figure 4E). The southernmost group consists of at least six flares rising from 1,560 m. In both groups, the acoustic signature of the rising bubbles reaches the depth of ∼750 mbsl. The total area of acoustic anomaly close to the seafloor measures 0.21 km². One flare from the southernmost group was used for bubble size and density parametrization, yielding gas flux estimates between 1.47 and 5.58 × 10⁶ kg/year for this site. The seismic data reveal a ∼400 m long, strong reflection with negative polarity right below the seabed at the location of gas seepage, indicative of widespread free gas in the shallow sediments (Figure 8C). Below the free gas accumulation, a column-shaped region of localised seismic blanking extends downwards towards the base of the GHSZ, in a region of disrupted reflections in the vicinity of an apparent BSR shoaling.

Urutī Ridge West

This is the shallowest of the analysed seep fields, and it lies in a region of relatively flat bathymetry at ∼1,140 m water depth. Urutī Ridge West is a SW-NE trending anticline that represents the southern extension of Urutī Ridge. The seismic profile shown in Figure 9 runs parallel to the strike of the anticline, and crosses two areas of gas seepage. The sedimentary sequence is characterised by relatively flat and parallel strata. The thickness of the GHSZ at Urutī Ridge West is ∼0.5 s, or ∼450 m using an estimated seismic velocity of 1800 m/s. While the BSR appears as a distinct negative polarity reflection adjacent to the seep locations, it is characterised by a series of lower amplitude reflections in the central part of the seismic profile, and it is not imaged beneath the regions of gas expulsion. High amplitude reflections with the same polarity as the seafloor probably point to the presence of concentrated gas hydrates or authigenic carbonates in the shallow sediments, while column-shaped regions of seismic blanking suggest upward fluid migration from the base of gas hydrate stability (BGHS) towards the seafloor. Similar to the Urutī Ridge eastern flares and to the Pahaua Bank northern flares, the gas flares at Urutī Ridge West are aligned roughly perpendicularly to the direction of plate convergence. The 18 imaged flares can be grouped into four clusters (Figure 4D), and they rise ∼500 m from the seafloor. One flare was selected to estimate bubble parameters (see Table 2). With a total area ∼0.45 km² of acoustic anomaly, the gas flux estimates for the Urutī Ridge West venting field lie between 1.82 and 6.54 × 10⁶ kg/year.

FIGURE 9. Overview of the Urutī Ridge West seep site: the interpreted TAN1808-02 seismic profile is shown in the main panel. The black arrows show the direction of fluid flow as interpreted from the seismic data and explained in the text. The black crossed circles represent the direction of fluid flow going into the page. BSR: bottom simulating reflection. (A) Map view of the Urutī Ridge West seeps (red box in Figure 4D), showing the acoustic backscatter anomaly in the echo integrated MBES (EM302) data in proximity of the seafloor. The location of the TAN 1808-02 seismic line and of the single-beam data (EK60) are indicated by the black and the red lines, respectively. The black dots show the locations of the main gas flares. The hydroacoustic data were collected during the TAN2012 voyage. (B) Echogram of a gas flare as imaged by the 38 kHz channel in the single-beam data. The horizontal scale is the same as the seismic panel.

Mungaroa Ridge

The Mungaroa seep field (Figure 10) is the deepest analysed in this study, with its main gas flare located at 2080 mbsl at the top of Mungaroa Ridge, a thrust-related four-way closure that lies at the toe of the accretionary wedge. Crutchley et al. (2021) investigated the gas hydrate system and fuid flow processes at Mungaroa Ridge in detail, using seismic reflection and multibeam data. They interpreted a gas-water contact pointing to a thick free gas column beneath the BGHS (Figure 10). This gas column is sufficiently thick to cause hydraulic fracturing through the gas hydrate stability zone, which is evidenced by a vertical chimney structure connecting the gas reservoir to the seafloor gas flare. Despite the existence of normal faults beneath the ridge, Crutchley et al. (2021) noted that they are not exploited for focused gas flow into the hydrate stability zone. We surveyed the region of gas seepage at the seafloor of Mungaroa Ridge during the three R/V Tangaroa voyages. From these data, we observed six flares rising from the seabed up to roughly 600 mbsl, making them the highest flares observed in the region (∼1,400 m high). The estimated methane fluxes at this site range from 2.52 to 7.18 × 10⁶ kg/year.

FIGURE 10. Overview of the Mungaroa Ridge seep site: the interpreted TAN 1808-97 seismic profile is shown in the main panel (interpretation after Crutchley et al., 2021). The bold black arrow shows the direction of fluid flow as interpreted from the seismic data and explained in the text. BSR: bottom simulating reflection; gwc: gas-water contact. (A) Map view of the Mungaroa Ridge seeps (red box in Figure 4A), showing the acoustic backscatter anomaly in the echo integrated MBES (EM302) data in proximity of the seafloor. The location of the TAN1808-97 seismic line and of the single-beam data (EK60) are indicated by the black and the red lines, respectively. The black dots represent the locations of the main gas flares. The hydroacoustic data were collected during the TAN2012 voyage. (B) Echogram of a gas flare as imaged by the 38 kHz channel in the single-beam data.

Discussion

The quantitative study of water column acoustic backscatter data combined with observations of subsurface geological structures has allowed a detailed characterisation of the five targeted cold seep areas on the southern Hikurangi Margin.

Sources of Uncertainty for Flux Estimations

There are several sources of uncertainty in the resulting estimations of gas fluxes at the seafloor.

The major source of uncertainty comes from the theoretical model used to predict bubble rising velocities. Because of the lack of quantitative observational data in the study area, we opted for using both clean and coated bubble models (Leifer and Patro, 2002) for our flux estimates.

Another type of uncertainty to be considered is related to instrumental parameters, that account for uncertainties of calibration, sound-velocity profiles and absorption rates. The calculated uncertainty for the SBES absolute backscatter measurements is ∼0.2 dB. The uncertainty of the absolute $S_v$ backscatter data will propagate on through the calculation of both BSD parameters and density of bubbles and will eventually affect flux estimations. However, due to the relatively high precision of the echo-sounders, we deem the effect of instrumental uncertainty negligible compared to the other sources of variability described in this and in the next section.

For the MBES data, usually 1–2% of the water depth is considered a conservative uncertainty in terms of positioning of soundings. Considering the deepest flare at Mungaroa Ridge (2080 mbsl) the highest uncertainty related to the spatial extent of the acoustic anomalies in the MBES data is of ∼1,600 m². This translates in ±2 × 10⁴ and ±3.01 × 10⁵ kg/year of methane considering coated and clean bubble models, respectively.

Uncertainties related to the estimations of bubble parameters, such as mean radius, bubble density and bubble size distribution, are addressed in detail in the next section.

Constraints on Bubble Size Distributions

Quantification of gas flux is dependent on observations of bubble parameters. Ideally, optical measurements such as video observation, bubble size measuring, and sampling of the seep fluids provide the most accurate measures of the bubble size distribution function, their rising velocity, and the chemical gas composition, enabling the determination of realistic values of gas flow rates (e.g., von Deimling et al., 2011; Higgs et al., 2019; Wang et al., 2016; Weber et al., 2014).

Due to the lack of optical observations of the seeps analysed in this study, no measurements of BSD are available, and we adopt a variation of the method proposed by Veloso et al. (2015) to estimate the BSD from the split-beam EK60 data. To test the validity of the results, we analyse the dependence of the estimated methane fluxes on different BSD functions: we first parametrise the BSD by assuming log-normal and Weibull probability density functions (PDF), and then compare the inverted results (Figure 11). The choice of these PDFs was made based on published seep studies, which have suggested several distribution functions to describe bubble size data including normal (Römer et al., 2012), log-normal (Veloso et al., 2015; Wang et al., 2016; Li et al., 2020), and Weibull (Dey and Kundu, 2012).

FIGURE 11. Bubble size distributions estimated from linear inversion of the split-beam data at the Mungaroa Ridge gas flare, imposing log-normal (A) and Weibull (C) distributions. Each curve is representative of a 50 m high horizontal slice of the gas flare, and is color coded according to the water depth. (B) and (D) show the theoretical frequency response curves at each horizontal slice of the gas flare, and the dots represent the observed S_v at the same water depth for the 18 and the 38 kHz channels.

The modelled $S_v$ depends on both the BSD and the density of bubbles (e.g., MacLennan et al., 2002). If the BSD parameters are over corrected during the curve matching procedure, the density of bubbles will be underestimated, so that the final impact on the calculation of gas fluxes is minimal. This can be observed, for example, in the inversion results for the Mungaroa Ridge flare at 2060 mbsl (Figure 10), where the mean bubble radius of the BSD modelled assuming a Weibull distribution increases above 2 mm, the density of bubbles for the corresponding depth shows a sharp decrease. The number of bubbles smaller than 0.5 mm is probably overestimated by the Weibull PDF, and the log-normal curves are more likely to represent the real BSD at these seep sites.

The use of two different approaches to estimate bubble parameters allows us to quantify the degree of agreement between the two methods. The fluxes estimated assuming a constant BSD for the entire flare and based on the normalised $S_v$ response (Approach 1–manual curve fitting method) show a general smooth trend with depth (Figures 5K–O). This trend can be explained by the fact that in this approach a mean $\overline{S_v}$ is calculated by averaging the echo integrated backscatter volume strength of the cells located at the same depth. Through this approach, we obtain one mean bubble radius per site rather than a profile of radii changing with depth. Because of the stronger response of bigger bubbles close to the emission site this method will also tend to overestimate the mean of the BSD, as the overall response of the flare is dominated by this area, as observed in Figures 5A–E. Despite this approximation, the manual curve fitting approach is a useful method that can be used to obtain informed starting models for the linear inversion method. Because of its deterministic nature, the second approach (linear inversion) requires the existence of a solid starting model in order to converge to a global minimum. This method provides a more reliable tool to track changes of physical parameters of the bubbles in the water column. The fluxes calculated through the linear inversion method (both assuming Weibull and log-normal distributions) show more variability in the calculated density of bubbles and fluxes, despite a general increase with depth (e.g., Figure 5). Table 2 summarizes the differences in the results of the two approaches. The differences between flux estimates calculated using a constant BSD for the entire flare (approach 1) and the fluxes obtained through the linear inversion process (approach 2) can be over 400%. For these reasons, we deem the flux estimates calculated through linear inversion more reliable (Table 2), and use those to extrapolate margin-wide flux estimates.

Source of Gas and Seismic Manifestation of Fluid Flow

The southern Hikurangi Margin is a well-established province of gas hydrate occurrence, focussed fluid migration and gas seepage (e.g., Barnes et al., 2010; Crutchley et al., 2019; Kroeger et al., 2019; Watson et al., 2020). Consistent with most subduction margins, the analysis of gases emitted at the seafloor suggests a predominantly microbial origin of methane over a thermogenic origin (Greinert et al., 2001; Faure et al., 2010). The co-location of high-resolution seismic reflection images and water column imaging (EK60 data) enables us to compare the sub-seafloor with the water column at gas seep locations.

At all the flare sites observed in this study, we found evidence of gas migrating from beneath the base of the GHSZ to the seafloor. In the case of Mungaroa Ridge, a large free gas reservoir in the core region of the anticline has been interpreted by Crutchley et al. (2021) as the source that supplies gas to the main seep observed at the seafloor. Here, Crutchley et al. (2021) suggest that over-pressured gas causes hydraulic fracturing of the overlying sediments leading to the formation of the vertical gas chimney imaged in the seismic data (Figure 10). The presence of such large and interconnected free gas accumulations is not observed at the other target sites of this study. However, highly reflective strata are imaged directly beneath the base of the GHSZ at Urutī and Glendhu ridges, as well as at Pahauau Bank. The enhanced reflectivity is likely to be caused by the strong impedance contrast between fine-grained low-permeability layers and sandy gas-charged sedimentary units, as interpreted by Turco et al. (2020) at Glendhu Ridge. Stratigraphically driven fluid migration along permeable dipping strata has been suggested to be the main mechanism for upward fluid flow in many anticline-related ridges on the Hikurangi Margin (e.g., Crutchley et al., 2019; Wang et al., 2017a; Turco et al., 2020; Barnes et al., 2010; Kroeger et al., 2021). Urutī Ridge is a good example of this process, where highly reflective strata appear to be transporting gas from depth upward into the GHSZ (Figure 6). At Urutī Ridge, Urutī Ridge West, Pahaua Bank ,and Glendhu Ridge, fluid migration through the GHSZ is also identified by areas of decreased amplitude (seismic blanking) and disrupted stratigraphic reflections beneath the seeps (Figures 6–9). Seismic blanking is often caused by strong signal attenuation at the seafloor or within the shallow sub-seafloor, caused by the presence of highly reflective interfaces. Such interfaces could come from authigenic carbonates or gas hydrate accumulations (e.g., Bohrmann et al., 1998). The disruption of reflections can be due to the scattering of seismic energy caused by the presence of gas (Judd and Hovland, 1992) or by physical disruption of layering caused by focused gas migration (e.g., Davis 1992; Gorman et al., 2002). While free-gas occurrence is likely to cause most of the seismic blanking observed beneath the seep sites, the presence of autigenic carbonates on the seafloor might contribute to the loss of seismic energy transmission in high-frequency data. In summary, the seismic images show a diversity of manifestations of free gas in the sub-seafloor, ranging from gas-water contacts at the base of a free gas reservoir (Mungaroa Ridge; Figure 10) through layer-parallel gas migration (e.g., Urutī Ridge; Figure 6) to vertical gas migration facilitated by hydraulic fracturing (Flemings et al., 2003; e.g., Pahaua Bank and Mungaroa Ridge; Figures 8, 10).

Temporal Variability of the Seeps

There are many mechanisms that control the activity of different types of gas seeps. Consequently, the time scales over which the activity of cold seeps fluctuates can span from minutes to millennia. For example, Feseker et al. (2014) document the eruption of a deep-sea mud volcano that triggered large methane and CO₂ emissions over a period of minutes. Pressure changes at the seafloor caused by tides have been shown to impact the flow rate of shallow and deep-sea gas seeps (Boles et al., 2001; Römer et al., 2016; Riedel et al., 2018), while seasonal sea-bottom temperature variations can cause cold seeps to hibernate during the cold months, trapping gas in the sediments that is released in pulses during warmer months (Berndt et al., 2014; Ferré et al., 2020). On the other hand, natural seismicity (Franek et al., 2017) and ocean warming (Baumberger et al., 2018) are potential triggers for significant release of methane from the sediments, especially in hydrate provinces.

Our study represents an opportunity to analyse the variability of methane emissions on the southern Hikurangi Margin over a three-year period. Although quantitative estimates were calculated only once for each gas flare, from either the 2019 or from the 2020 datasets, from the qualitative analysis of multibeam and split-beam data, no substantial difference could be observed in the activity of the main seeps at the time of each survey (Figure 12). In fact, the spatial extent of the acoustic anomaly close to the seafloor remains constant for the five target areas in the three datasets, as does the height of the acoustic flares.

FIGURE 12. Evolution of the Glendhu Ridge seep site over the years. The panels on the left (A,C,E) show range stacked views of MBES (EM302) data of the gas flares on the top of Glendhu Ridge from TAN1808, TAN1904, and TAN2012 datasets, respectively. The panels on the right (B,D,F) show the acoustic backscatter anomaly in the echo integrated MBES data in proximity of the seafloor. The black lines represent the ship track shown in the left panels. The coordinate system is UTM Zone 60 S (WGS84 datum).

In addition to acoustic observations, it is known from authigenic carbonates (for example on Urutī Ridge) that many of the seep sites have been active for thousands of years (e.g., Jones et al., 2010; Liebetrau et al., 2010). Likewise, there is evidence for stable methane seepage over intermediate timescales from tube worms (Lamellibrachia spp.) up to 2 m long sampled at Mungaroa, Urutī and Glendhu ridges (TAN1904 Voyage Report, NIWA). Tube worms of this species require at least 200 years to reach such lengths (Fisher et al., 1997; Cordes et al., 2007).

While the combination of acoustic observation and bio-geological sampling might indicate a constant seepage activity throughout this time, we cannot rule out that methane fluxes vary over seasonal or shorter cycles. For simplitity, we assume a constant discharge rate for the flux estimates presented in this work. Understanding and monitoring the temporal variability of a field of cold seeps is relevant to several scientific and socio-economic issues. At a national scale, one of the most direct implications is related to regional ecosystem management. Cold seeps are increasingly recognized as centres of local biogeochemical cycling and oases for many animals with recent studies finding that commercially important fisheries species are associated with seep habitats and consume methane derived carbon from chemosynthetic production in seep systems (Grupe et al., 2015; Levin et al., 2016; Seabrook et al., 2019; Turner et al., 2019).

Margin-Wide Estimates of Seafloor Methane Flux

To refine our understanding of the global carbon budget, it is important to study the potential implications of seabed gas release on a regional scale. To this end, the relevance of margin-wide studies on natural methane seeps has increased in the past decade: Pohlman et al. (2011) found that up to 28% of the total dissolved organic carbon derives from fossil methane, while Garcia-Tigreros et al. (2021) conclude that aerobic oxidation of CH₄ has a greater influence on ocean chemistry in regions where methane concentrations are locally elevated. Based on the analysis of more than 300 gas seeps, Riedel et al. (2018) estimate a combined average in-situ flow rate of about 88 × 10⁶ kg/year for the Cascadia Margin. Sahling et al. (2014) calculate the gas flux on the western margin of Svalbard to be between 0.11 and 1.89 × 10⁶ kg/year, while the combined methane leakage from the seafloor on the US Atlantic Margin is estimated to range between 1.5 and 90 × 10⁴ kg/year (Skarke et al., 2014). About 0.13–1.01 × 10⁶ kg of methane per year have been estimated to leak from gas seep fields on the Makran Margin offshore Pakistan (Römer et al., 2012).

The distribution of gas seeps on the Hikurangi Margin has been investigated in detail by Watson et al. (2020), who identified 1,457 gas flares from water column data, spanning from East Cape to Kekerengu Bank, off Kaikōura peninsula, in an area of approximately 51,215 km² (Figure 2A). If we consider the flux estimates presented in this paper for the five analyzed seep fields, and we divide them by the number of flares observed at each site, we obtain average methane fluxes per flare between 1.9 × 10⁵ and 6.4 × 10⁵ kg/year, considering coated bubbles and clean bubbles models, respectively. Multiplying by the number of observed flares from Watson et al. (2020), we can extrapolate a total flux ranging between 2.77 × 10⁸ and 9.32 × 10⁸ kg of methane released each year on the whole Hikurangi Margin. These estimates correspond to, respectively, ∼67 and ∼233% of the total amount of methane that was released in 2019 from land sources in New Zealand (i.e., 4.121 × 10⁸ kg/year New Zealand Ministry for Environment—Manatū Mō Te Taiao, www.environment.govt.nz). It is worth pointing out that these estimates might represent an underestimation of the real methane flux occurring on the margin, because they are based solely on acoustic evidence of focused fluid flow (i.e., gas bubbles in the water), and do not take into account diffusive fluid seepage, which is also inferred to occur on the Hikurangi Margin (Watson et al., 2020). On the other hand, assuming a constant methane discharge at the seep sites over seasonal or shorter cycles might lead to an overestimation of annual gas fluxes. Similarly, the assumption that all flares in the same seep field have identical bubble parameters is an approximation that can result in overestimating methane discharge at smaller seeps. To better assess the methane fluxes, long term observational data and chemical sampling at the analysed sites might be required.

Potential Implications for Ocean pH and Local Deoxygenation

From the acoustic imaging of the water column, no evidence was found that methane bubbles reach the sea surface at any of the analysed flares. Considering the depth of the seep sites reported here, two different processes will interact to prevent the CH₄ emitted at the seafloor from reaching the atmosphere: 1) the methane contained in the bubbles will dissolve into the water driven by the low concentration of CH₄ in the ocean (Wiesenburg and Guinasso, 1979), and2) dissolved CH₄ is converted into CO₂ in the water column by abiotic and biotic forces (McGinnis et al., 2006). Although bubble-stripping and methane oxidation reduce the amount of CH₄ released into the atmosphere, these processes can significantly impact local marine habitats and ocean chemistry.

As the majority of the CH₄ emitted by cold seeps remains in deep waters, aerobic oxidation is a primary sink for the methane, as well as a source of CO₂. This CO₂ production needs to be considered with respect to acidification of deep water (Archer et al., 2009; Biastoch et al., 2011), particularly as seeps represent a more direct source in deep water than the transfer of anthropogenic carbon via deep water formation and transport. Large-scale methane release has resulted in ocean acidification in earths geological past (Zachos et al., 2005), but assessments of the current contribution of methane seeps indicate a relatively minor impact on deep water pH (Garcia-Tigreros and Kessler, 2018). However, this source may become significant in response to warming and associated destabilisation of methane hydrates in regions such as the Arctic (Biastoch et al., 2011).

In the current study we estimated the regional contribution to deep waters on the Hikurangi Margin by scaling up the methane release estimated from the bubble plumes to the total number of flares identified in the acoustic data (Watson et al., 2020). Assuming this methane loading was uniformly distributed within the bottom 100 m of water overlying the sediment for the entire margin, from East Cape to the Kaikoura peninsula, provided an estimate of total input of 2.77 × 10⁸ and 9.32 × 108 kg methane/year, based upon the clean and coated bubbles models, respectively. If 100% of this methane is oxidised to CO₂ then the resulting change in the carbonate system, calculated using measured bottom water dissolved inorganic carbon, total alkalinity, temperature, and salinity (C. Law, pers. comm.) in the CO₂sys programme (Hunter, 2015), would result in a decrease in pH_T of 0.048–0.144 relative to the background value of 7.962. These estimates are conservative, and suggest a relatively minor impact to a significant decrease in regional pH, with the upper estimate exceeding the surface ocean pH decrease arising from anthropogenic CO₂ emissions to date (Orr et al., 2005), and comparable to 50% of the projected pH decrease in New Zealand waters by the end of this century (Law et al., 2018). A reduction in pH of this magnitude could have significant impacts on benthic calcifying organisms, particularly as the Hikurangi Margin gas seeps are within the regional depth range of the Aragonite Saturation Horizon (Bostock et al., 2015), below which solid carbonate becomes thermodynamically unstable. However, this estimate of pH decrease should be regarded as an upper limit, as dilution and transport of methane in and out of the region are not considered in this estimate. Refined estimates of pH change require direct measurement of dissolved methane and regional modelling of methane distribution and dispersion using ROMS (Hadfield et al., 2007).

The aerobic oxidation of methane (methanotrophy) by free-living and symbiont-associated bacteria occurs in the seafloor and water column around sites of gas release (Steinle et al., 2015; Sweetman et al., 2017; Levin, 2018). Although not analysed in this study, the aerobic oxidation of methane to CO₂ can also cause localized deoxygenation at a regional scale, impacting ecosystem health and species distribution (Boetius and Wenzhöfer, 2013; Breitburg et al., 2018). Throughout the global ocean there has been an oxygen loss of at least 2% over the past 50–100 years, with well understood linkages of ocean warming reducing oxygen solubility and ocean ventilation (Dickens, 2001; Levin, 2018). However, the influence of features that influence dissolve oxygen at regional scales, such as methane seeps, remains unclear. In addition to the direct effect of methanotrophy, some evidence suggests that seeps with strong bubble plumes, such as observed at some of sites reported here, can draw nutrient and hydrocarbon rich water towards the surface, stimulating primary production and eventually drawing down oxygen as well (Levin, 2018). The potential scale of this process could be significant, as seen in the aftermath of the Deepwater Horizon oil spill in the Gulf of Mexico, where a major reduction of oxygen was measured ad microbial communities released from the spill (Kessler et al., 2011). The downstream ecological and biogeochemical impacts of the release of large volumes of methane by the Hikurangi Margin seeps, as detailed in this study, warrants further scientific attention to understand the implications of future change.

Conclusion

The combination of seismic and hydroacoustic data analysis allowed the characterisation of five cold seep sites on the southern Hikurangi Margin in terms of geological setting and gas flux estimates. Seismic imaging of the geological structures underlying the seep sitesprovided insights into the origin of the gas in the subsurface. Hydroacoustic data collected over 3 years allowed mapping of the backscatter anomalies near the seafloor at the sites of seepage and pinpointing the location of the main gas flares on the seabed. A total of 43, 33, and 53 individual flares were identified from the TAN1808, TAN1904, and TAN2012 datasets, respectively.

The use of the multi-frequency split-beam echosounder allowed estimates of the gas flux rates at the five target sites to be made. The five cold seep fields analysed in this study on the southern Hikurangi Margin of New Zealand lie in water depths ranging from 1,110 to 2060 m, and emit, combined, between 8.66 and 27. 21 × 10⁶ kg of gas per year. The extrapolated methane flux for the whole Hikurangi Margin range between 2.77 × 10⁸ and 9.32 × 10⁸ kg of methane released each year. These estimates are based on acoustic evidence of focused fluid flow, and do not take into account diffusive seafloor seepage.

The results of this study provide the most quantitative assessment to date of total methane release on the Hikurangi Margin, filling gaps of unknown methane sources and better constraining models of ocean acidification and deoxygenation. Moreover, the fluxes presented here can be used as a proxy to monitor changes in the flux rates over the mid-to long term associated with ocean warming.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

FT: conceptualization, literature review, acoustic, and seismic data acquisition (TAN1808–TAN2012), acoustic data processing, software developing for flux estimates, writing. YL: gas bubble project leader, supervision in the methodology, software developing for acoustic data analysis, and processing, data curation. SW: gas bubble project leader, writing–review, supervision, seismic, and acoustic data acquisition (TAN1808). CL: took care of the ocean acidification implications and estimates on pH change, writing–review. GJC: gas hydrate seismic project leader, conceptualization of manuscript, especially discussion, seismic and acoustic data acquisition and processing (TAN1808), writing–review. SS: contribution to writing biogeochemical implication of the study, acoustic data acquisition (TAN1904). JM: supervision, seismic and acoustic data acquisition (TAN1808), supervision, project management and facilities. CL: contribution to calculating CO₂ for ocean acidification, acquisition of acoustic data (TAN1904). IP: gas bubble project co-leader, acoustic data acquisition (TAN2012), writing–review. JITH: gas hydrate project leader, seismic and acoustic data acquisition (TAN1808), writing–review. SW: seismic data acquisition and processing (TAN1808), acoustic data acquisition (TAN2012). AG: supervision, provided facilities and funding.

Funding

Part of this work has been carried out within the HYDEE research programme, funded by the New Zealand’s Ministry for Business, Innovation and Employment, contract no. C05X1708, and is included a Ph.D. thesis from the University of Otago (Turco, 2021). Part of the funding came from the Smart Idea project Broadband acoustic characterization of free gases in the ocean water, contract nr. C01X1915.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References