Multiple drivers and controls of pockmark formation across the Canterbury Margin, New Zealand

Shallow seabed depressions attributed to focused fluid seepage, known as pockmarks, have been documented in all continental margins. In this study, we demonstrate how pockmark formation can be the result of a combination of multiple factors—fluid type, overpressures, seafloor sediment type, stratigraphy and bottom currents. We integrate multibeam echosounder and seismic reflection data, sediment cores and pore water samples, with numerical models of groundwater and gas hydrates, from the Canterbury Margin (off New Zealand). More than 6800 surface pockmarks, reaching densities of 100 per km2, and an undefined number of buried pockmarks, are identified in the middle to outer shelf and lower continental slope. Fluid conduits across the shelf and slope include shallow to deep chimneys/pipes. Methane with a biogenic and/or thermogenic origin is the main fluid forming flow and escape features, although saline and freshened groundwaters may also be seeping across the slope. The main drivers of fluid flow and seepage are overpressure across the slope generated by sediment loading and thin sediment overburden above the overpressured interval in the outer shelf. Other processes (e.g. methane generation and flow, a reduction in hydrostatic pressure due to sea‐level lowering) may also account for fluid flow and seepage features, particularly across the shelf. Pockmark occurrence coincides with muddy sediments at the seafloor, whereas their planform is elongated by bottom currents.


| INTRODUCTION
Fluid seepage can significantly shape the seafloor across continental margins, giving rise to a range of morphologies. Pockmarks-seabed depressions with circular to elongate planforms, steep flanks and flat to cone-shaped bottoms formed by focused fluid flow and escape-are amongst the most ubiquitous of these morphologies (Andresen & Huuse, 2011;Fader, 1991;Gay & Berndt, 2007;Judd & Hovland, 2007;Solheim & Elverhøi, 1985). Pockmarks provide an insight into sub-seafloor plumbing systems and basin dynamics (Talukder, 2012), and are thus useful in hydrocarbon exploration and carbon sequestration assessment (Hovland & Judd, 1988). Fluid emissions at pockmarks are a driver of cold seep ecosystems (Foucher et al., 2015;Levy & Lee, 1988) and may be linked to past and present climate change Westbrook et al., 2009). Since pockmarks are often found in the vicinity of fluid-driven sedimentary failures, their study is also important for hazard assessment (Deville et al., 2020;Hovland et al., 2002;Sills & Wheeler, 1992).
Depressions located off the east coast of the South Island of New Zealand have been the focus of intensive research in recent years. Most of the attention has been paid to the Chatham Rise (Figure 1a), which hosts depressions ranging 0.15-12 km in diameter in a depth range of 500-1100 m. Davy et al. (2010) had associated pockmarks located in the Chatham Rise, in 500-700 m water depth and with diameters of up to 150 m, with gas hydrate dissociation triggered by sea-level fall during glacial periods and higher bottom-water temperatures during interglacial warming. Their hypothesis was based on the coincidence of the upper water depth limit of pockmarks on the seafloor with the upper limit of the methane hydrate stability zone in the ocean. However, this inference was later refuted due to the absence of both methane in sediment pore water and bottom simulating reflections (BSR) in seismic reflection data (Bialas et al., 2013;Klaucke et al., 2018). Klaucke et al. (2018) and Waghorn et al. (2018) investigated depressions of the Chatham Rise located at water depths of 600-1100 m and with diameters of up to 10 km. They reported widespread polygonal faults that may have acted as fluid migration pathways in the past. The depressions are thought to have formed by scouring and deposition by strong hydrostatic pressure due to sea-level lowering) may also account for fluid flow and seepage features, particularly across the shelf. Pockmark occurrence coincides with muddy sediments at the seafloor, whereas their planform is elongated by bottom currents.

K E Y W O R D S
Canterbury Margin, groundwater, methane, pockmark, sediment loading Highlights • We document >6800 pockmarks across the Canterbury Margin, with densities reaching 100 per km 2 • Pockmarks were formed by expulsion of methane, and saline and freshened groundwater • Overpressure due to sediment loading is the main driver of fluid flow across the slope • Methane generation and lower sea-levels may account for fluid flow across shelf • Pockmark distribution and shape are influenced by muddy seafloor sediments and bottom currents bottom currents at nucleation points generated by fluid venting. Klaucke et al. (2018) suggested that the fluid involved originated from sediment dewatering during opal A/CT transformation. More recently, Stott et al. (2019) suggested that depressions across the Chatham Rise formed due to release of 14 C-dead carbon dioxide (CO 2 ) and carbon-rich fluids from subsurface reservoirs, with the most likely source being dissociated Mesozoic carbonates that subducted during the Late Cretaceous. Depressions near the Otago submarine canyon complex (Figure 1a), on the other hand, were investigated by Hillman et al. (2015). Here, the largest concentration of depressions is reported adjacent to actively eroding canyons, and their formation is attributed to groundwater flow and interaction with the Southland Current.
The importance of bottom current systems in reshaping F I G U R E 1 (a) Bathymetric model of the eastern continental margin of the South Island of New Zealand (Source: GEBCO). The inset shows the map of New Zealand and the location of (b). (b) Topography of the Canterbury Plains and offshore margin (Source: National Institute of Water and Atmospheric Research (NIWA)). The location of the Sub-tropical Front (STF) is shown . The transparent grey zone shows seafloor covered with sediments containing >50% mud (Bostock et al., 2019) depressions in both the Chatham Rise and Otago Margin is highlighted by Hillman et al. (2018).
In this study, we focus on the Canterbury Margin (Figure 1b), which is located between the Chatham Rise and Otago Margin, and which is the least explored sector of the eastern passive margin of the South Island in terms of fluid flow and seepage. By integrating geophysical and core data with numerical modelling, we: (i) document widespread evidence for fluid flow and seepage, (ii) identify the types of fluids involved and (iii) infer the drivers and controls of pockmark formation. We demonstrate how pockmark occurrence, characteristics and formation can be the consequence of the flow of different fluid types, overpressures, seafloor sediment type, stratigraphy and the influence of bottom currents.

| REGIONAL SETTING
The 50,000 km 2 Canterbury Basin is part of the eastern continental margin of the South Island of New Zealand, which includes the Canterbury braid plains onshore and the Canterbury Bight and continental slope offshore ( Figure 1b) (Browne & Naish, 2003). The continental shelf extends about 180 km from south-west to north-east and reaches ca. 95 km from the coastline, with a maximum gradient of 0.1°. Beyond the shelf break (ca. 140 m below sea-level (bsl)), the gradient increases to 2° and >5°, for the north-eastern and south-western sections of the slope respectively. Consequently, the width of the continental slope changes from 30 km in the north-east to 14 km in the south-west. The base of the slope is located at a depth of 1100 m bsl (Browne & Naish, 2003;Lu & Fulthorpe, 2004).
The Canterbury Basin forms part of a passive margin located on the landward edge of a continental fragment that separated from Marie Byrd Land (West Antarctica) ca. 100 million years ago (Ma), after the break-up of the Gondwana supercontinent (Strogen et al., 2017). Since then, tectonic activity offshore has mainly been associated with subsidence in the basin centre, resulting in limited fault activity during the Cenozoic (Brown & Field, 1988;Lu et al., 2003). Convergence of the Pacific and Australian plates at the end of the Oligocene (ca. 23 Ma) led to uplift of the Southern Alps, leading to enhanced erosion and high sediment supply to the Canterbury Basin (Lu et al., 2003).
The sedimentary architecture of the Canterbury Margin was formed under the influence of a long-term relative sea-level cycle starting ca. 80 Ma (Fulthorpe et al., 1996;Lu & Fulthorpe, 2004;Lu et al., 2005). Lu and Fulthorpe (2004) identified at least eight complete depositional sequences (S12 to S19) in the stratigraphic record of the last ca. 2.5 My. Some sequences, dated to the early Quaternary, are characterised by seismic morphologies indicative of sediment drifts associated with palaeo-slopes (Lu et al., 2003). Sedimentation rates during the Quaternary average ~25 mm a −1 , although values of ca. 150 mm a −1 were reached 0.25 Ma (Lu et al., 2005;Villasenor et al., 2015). Lithological descriptions from International Ocean Discovery Program (IODP) Expedition 317 cores show a dominance of mud and sandy mud interrupted by sand layers (Fulthorpe et al., 2011). The surficial sediment distribution across the margin has been mapped by Bostock et al. (2019). The inner shelf, especially in proximity to the mouths of the Ashburton and Rakaia Rivers, predominantly consists of sand and gravel. The middle to outer shelf is covered by sand in the north-eastern section and mud in the south-western section (Figure 1b). The slope is predominantly covered by mud.
Two surface water currents flow off the east coast of the South Island: the southward-flowing sub-tropical East Cape Current and the north and eastward-flowing extension of the Southland Current, which carries sub-Antarctic waters (Chiswell, 2003;Sutton, 2001). These currents are separated by the Sub-Tropical Front (STF), a complex and irregular zone of enhanced meridional temperature and salinity gradients that encircles the globe at approximately 45°S (Heath, 1985;Sutton, 2003) ( Figure 1b). Surface waters of the STF flow north-east along the continental margin before being deflected by the Chatham Rise at a latitude of 43°S; the strongest thermal gradients occur at approximately 500-m water depth (Chiswell, 1994;Chiswell & Rickard, 2006;Sutton, 2003). Water depths between 600 m and 1450 m are dominated by the northward-spreading Antarctic Intermediate Water, whereas the Circumpolar Deep Water flows below 1450 m McCave & Carter, 1997).

| Multibeam echosounder (MBES) data
MBES data were acquired using Kongsberg EM2040 and EM302 systems with frequencies of 200-400 kHz and 30 kHz respectively (Figure 2b). The systems were combined using a single Applanix POSMV V5 unit for online roll, pitch, heave and yaw corrections. Sound velocity profiles were acquired at regular intervals and employed during post-processing for calibration of travel-times and water depth. Two bathymetric grids, with a cell size of 2.5 m × 2.5 m on the shelf and 20 m × 20 m on the slope, were generated. In addition to bathymetry, water column backscatter data were acquired during TAN1703. These data were analysed visually for water column anomalies using QPS Fledermaus Midwater software.
On the shelf, the backscatter intensity in the multibeam data was extracted and analysed using the QPS Fledermaus software packages. Attenuation coefficients of 50 dB km −1 and 70 dB km −1 for 200 kHz and 300 kHz, respectively, were derived according to Ainslie and McColm (1998). The raw backscatter time series data for each footprint (snippets) were used, where available, and compensated for the angle of incidence with a reference grid. The final mosaic was created with a 1 m cell size.

| TOPAS sub-bottom profiles
During the TAN1608 and TAN1703 surveys, a TOPAS PS 18 Parametric Sub-Bottom Profiler (SBP) was deployed ( Figure 2b). Linear Chirp waves of 2.0-6.0 kHz frequency and a Chirp length of 15 ms/20 ms were used to acquire sub-bottom profiles with a vertical resolution of up to 20 cm and a penetration of up to 80 m. A motion data feed from the POSMV allowed stabilisation of the beam for heave, roll and pitch. In areas with steep seafloor gradients, the acoustic beam was steered manually.

| Multichannel seismic (MCS)
reflection data TAN1703 acquired 600 km of high-resolution MCS reflection data using a mini GI-gun (13/35 cubic inch), deployed at 1.5 m water depth, and shooting with a pressure of 1800 to 2000 PSI (124 to 138 bar) (Figure 2a). The receivers were within a 300 m active solid-state GeoEel digital seismic streamer (Geometrics), with a hydrophone group spacing of 12.5 m. The acquisition parameters were set to a shot interval of 3 s and a recording length of 1.5 s at a navigation speed of 4 knots, resulting in a shot spacing of 6 m. In general, the quality of the data was good, although noise levels were occasionally high due to unfavourable weather conditions.
Data processing was carried out using Globe Claritas and included the following operations: conversion from SEGD to SEGY, coordinate conversion, definition of streamer geometry, band pass filtering (corner frequencies of 50, 100, 500 and 700 Hz) and spherical divergence F I G U R E 2 (a) Spatial coverage of multichannel seismic (MCS) reflection data (TAN1703 and EW00-01) and location of TAN1703 sediment cores and scientific ocean drilling sites (ODP Leg 181, IODP Expedition 317). (b) Spatial coverage of TAN1608 and TAN1703 multibeam echosounder (MBES) data and sub-bottom profiles (SBP). Isobath interval is 50 m corrections, common depth point (CDP) binning (6.25 m bin size) and sorting, normal moveout correction, stacking and quality control of processed data. Post-processing included swell correction, deconvolution, migration and automatic gain control and horizontal trace balancing. The processed seismic profiles have a vertical resolution of 2-2.5 m.
All seismic data presented in this study conform with the American polarity convention, where an increase in acoustic impedance (e.g. at the seafloor) leads to a positive polarity peak in the seismic data (Brown & Abriel, 2014). The polarity of reflections beneath the seafloor can be difficult to discern, due to tuning of the wavelet caused by 'thin' layers (i.e. layers much thinner than the seismic wavelength). Impedance increases and impedance decreases are interpreted based on whether the strongest part of the wavelet at a given reflection is a peak (positive polarity) or a trough (negative polarity) respectively.

| Sediment cores
The seafloor was sampled using a purpose-built piston coring system during TAN1703. Six cores were recovered, two from the shelf (8, 26) and four from the slope (28-31) ( Figure 2a; Table S1). An additional 25 attempts were made at sampling the shelf, with no recovery, which is likely due to the occurrence of coarse sediments at the seafloor. The cores were cut into 1 m sections for storage, logging and geochemical analysis of sediment and pore water.
3.1.5 | Other data MCS reflection data collected during RV Maurice Ewing cruise EW00-01 (Lu et al., 2003) were utilised. These MCS data were acquired in 2000 and cover an area of ca. 4800 km 2 along the outer shelf and slope ( Figure 2a). The survey includes 57 profiles (ca. 3750 km in total) with a spacing of 0.3-0.7 km between the individual lines. For data acquisition, two GI air guns (45/45 in 3 ) and a streamer containing 96-120 channels in 12.5 m groups, each containing 26 hydrophones, were deployed. A shot interval of 5 s and a record length of 3 s were set during acquisition.
Bathymetric data with 100 m grid size for the entire Canterbury Margin were made available by NIWA, whereas borehole data were obtained from IODP Expedition 317 (Fulthorpe et al., 2011) and Ocean Drilling Program (ODP) Expedition 181 (Carter et al., 1999) (Figure 2a).

| Mapping of seafloor depressions
The ESRI ArcMap semi-automatic tool by Gafeira et al. (2018) was employed to map seafloor depressions in the MBES data. The output shapefile includes location, size, depth and shape of the depressions.

| Sediment core analyses
The analysis of TAN1703 sediment cores included visual description and elemental analysis via X-ray fluorescence (XRF) using a CORTEX core scanner. For interpretation we focused on three element ratios: Mn/Al, Ca/Ti and Sr/ Ca. The Mn/Al ratio serves as a proxy for redox conditions in marine sediments (Jaccard et al., 2009). Ca/Ti and Sr/ Ca are used for the identification of authigenic carbonate formation in the presence of methane (Li et al., 2016). Ca/ Ti ratios also provide a proxy for sediment origin (marine or terrigenous) (Croudace & Rothwell, 2015).
The measurement of anion and cation concentrations was carried out by Hill Laboratories in Hamilton, New Zealand, via Inductively Coupled Plasma Mass Spectrometry (IPC-MS) based on Rice et al. (2012). Samples for Cl − anions were prepared by ferric thiocyanate colorimetry (standard method APHA 4500 Cl − E). SO 4 2− anions were analysed via ion chromatography according to APHA 4110 B standard. Stable isotope measurements of δ 18 O, δ 2 H and δ 13 C DIC were carried out using a GasBench II coupled to a DELTA V Plus isotope ratio mass spectrometer (IRMS) (Thermo Fisher Scientific, Bremen, Germany) at NIWA. Sample treatment for δ 18 O and δ 2 H measurements followed methods described in Duhr and Hilkert (2004). Values of δ 13 C DIC for water samples were calculated from δ 13 C (g) values by taking into account the carbon isotope fractionation associated with the CO 2(g) -CO 2(aq) partition, as proposed by Assayag et al. (2006). Measured values were normalised using National Institute of Technology and Standards (NIST) reference materials. For absolute DIC concentrations in the pore water, coulometric titration (UIC CM5017) was performed in the isotope core laboratory of Texas A&M University. Resulting values were normalised with a multipoint calibration using in-house N 2 CO 3 solutions.

Groundwater flow
The governing equation describing groundwater flow used in this study is given as follows: where S s is the specific storage, is the sea-level elevation, k is the permeability tensor, µ is the fluid viscosity, g is the gravity constant, h is the freshwater hydraulic head, ω is the sedimentation rate, f is the fluid density and s is the sediment density.
Equation (1) was solved using the finite element method using the basin simulator RIFT2D (Person & Garven, 1994). Triangular elements using linear shape functions were employed. The resulting system of algebraic equations was solved directly using Gaussian elimination. The accuracy of using Equation (1) for representing compaction driven flow was tested by comparing the RIFT2D output to the analytical solution of Gibson (1958). RIFT2D results were within ca. 3% of the analytical solution.
The model domain is 173 km in length and 1.1-1.9 km in depth ( Figure S1). The subsidence rate varied between 0.23 and 0.55 mm a −1 (Browne & Naish, 2003), generating sediments 161-350 m thick over the 700 ka simulation period. Five lithofacies were included in our model ( Figure  S1). The distribution of the facies was loosely based on stratigraphy presented by Browne and Naish (2003). Table  1 lists the properties assigned to the five hydrostratigraphic units. Values of fluid viscosity, fluid density and sediment density were kept constant (Table 1). Specific storagethe aquifer's capacity to release water from or take into storage when water level changes-was varied between 3 × 10 −5 and 5 × 10 −6 m −1 , depending on lithology. The values of specific storage listed in Table 1 are typical for sand/gravel and conservatively low for clays (Konikow & Neuzil, 2007). As a result, our head value calculations, described below, are conservative. Porosity was allowed to decay with depth using a modified Athy's law expression based on effective stress and sediment compressibility: where ϕ o is the porosity at the sediment-water interface, σ e is effective stress and β is sediment compressibility.
In all simulations, sea-level was varied along the top boundary, modifying the imposed hydraulic heads for nodes below sea-level. No-flux boundary conditions were imposed along the sides and base of the model domain. For all nodes below sea-level, the following boundary conditions were imposed: where h bc (x, t) is the specified freshwater hydraulic head at the top of each nodal column, d (x, t) is the water depth below sea-level at a given location and time for a given sediment column, (t) is sea-level elevation and o is the base density at standard temperature, concentration and pressure (10°C, 0 mg/L, atmospheric pressure). Land surface nodes above sealevel were assigned a head value equal to the surface elevation.
Four scenarios were simulated (Table 2) to assess the influence of permeability anisotropy, sediment loading and sea-level loading on computed heads and pressures. Because of the lenticular nature of sediment grains and the presence of discontinuous fine-grained horizons within stratigraphic layers, groundwater tends to flow parallel to bedding. Permeability anisotropy accounts for this preferential flow parallel to bedding planes. Permeability anisotropy (k x /k z ) of the sediments (where k x and k z are the permeability parallel and perpendicular to bedding respectively) was varied between 30 and 50 to reproduce 'high' and 'low' permeability scenarios. These values are in the same order of magnitude as those used by Micallef et al. (2020) for the Canterbury Bight, and are reasonable for metre-scale clastic deposits (Sanford et al., 2004). In these simulations, both sedimentation and sea-level fluctuations were imposed as mechanical loads. In a third simulation, only sediment loading was imposed using 'low permeability' conditions. In the fourth scenario, only sea-level loading in 'low permeability' conditions was considered. The rates of sea-level change are up to an order of magnitude larger than sedimentation rates ( Figure S2). However, the average sea-level was close to zero over seven glacial cycles. The response time: Abbreviations: S s , specific storage; ϕ o , porosity at the sediment-water interface; β, sediment compressibility.
for the clay facies in our model was estimated at ca. 1 Ma for the low permeability scenario (where H, the sediment thickness, was assumed to be 1000 m). This is the time required for hydraulic heads to equilibrate to sediment and sea-level loading conditions in low permeability environments. Because of the long response time, computed heads were initialised assuming steady-state conditions. The model was then run for 700 ka using time steps of 1 ka.

Gas hydrate stability
The CSMHYD code from the Colorado School of Mines was used to model the gas hydrate stability zone in sediments on the slope (Sloan, 1998). CSMHYD predicts the thermodynamic stability of gas hydrate for a given pressure, temperature and gas composition. We assumed two different gas compositions to predict gas hydrate stability: (i) 100% methane (CH 4 ) and (ii) 100% carbon dioxide (CO 2 ). The geothermal gradient was estimated at 46°C km −1 at IODP Expedition 317 Site U1352 (Fulthorpe et al., 2011). However, in view of uncertainties during temperature measurements due to an irregular cooling curve, this value is considered unreliable. We therefore tested a range of 36-46°C km −1 . A pore water salinity of 3.5 weight % was used, and the bottom water temperature profile was obtained from Conductivity Temperature Depth (CTD) data available from the world ocean database (2018) and hosted by the National Oceanographic Data Centre (NODC) (https://data.nodc.noaa.gov).
Our approach was as follows: 1. Convert the seafloor depth grid (bathymetry) to a seafloor temperature grid using the NODC temperature profile ( Figure S3). 2. Calculate a seafloor pressure grid assuming a water density of 1027 kg/m 3 . 3. For each pressure and temperature value on the seafloor grid, calculate whether gas hydrate is stable under the assumption of (i) 100% methane as the hydrate-forming gas and (ii) 100% carbon dioxide as the hydrate-forming gas. From these grids, we determine two seafloor contours: (i) a contour that marks the upper limit of the pure methane hydrate stability zone within the sediments; (ii) a contour that marks the upper limit of the pure carbon dioxide hydrate stability zone within the sediments.

| Bathymetry
We have mapped >6800 depressions across the study area ( Figure 3a). Of these, ca. 2800 are located in the middle to outer shelf (80-140 m bsl) and have a density of up to 100 depressions per km², with the peak in depression density occurring at a depth of 105 m (Figure 3b). Estimates of density are based on the assumption that the multibeam swaths are characteristic of the area in general and depression distribution in the areas between the swaths is similar to that within them. Across the shelf, most of the features are shallow (<1 m deep) and have an average diameter of ca. 20 m. The shape of the depressions varies from circular to highly elongate. In the case of elongate depressions, the long axis is generally oriented parallel to the coastline (south-west to north-east), with a gentle slope on the south-western side and a steeper slope on the north-eastern side (Figures 3a and 4a). Circular and elongate depressions are evenly distributed across the shelf, although elongate depressions are somewhat more common in the vicinity of the shelf break.
Approximately 4000 depressions are located on the lower continental slope (500-1100 m bsl), reaching a maximum density of 35 depressions per km² (Figure 3a). The peak in depression density on the lower continental slope occurs at ca. 800 m bsl (Figure 3b). The depressions are circular to elliptical in planform, with a diameter of 50-200 m and a depth of up to 20 m. All depressions on the slope have a west-south-west to east-north-east orientation and an asymmetry that is opposite to that of the shelf depressions (Figures 3a and 4c). Very few depressions occur on the upper continental slope at depths of 140-500 m bsl.

| Backscatter
The depressions show variable backscatter strengths. While most of the depressions have no distinct backscatter difference compared to their surroundings, ca. 5% of depressions on the shelf show higher backscatter strength (ca. −15 dB difference compared to their surroundings). The elevated backscatter anomalies are predominantly associated with circular depressions, whereas the most elongate depressions tend to have similar backscatter strength to the surrounding seafloor ( Figure 4).

| Water column
We are able to identify numerous anomalies in the water column data. Less than 10 of these anomalies share the characteristics of gas flares. They are vertical to slanting, display strong backscattering with a high height to width ratio, and are connected to the seafloor. Such flares produce backscattering of the transmitted pressure wave due to the impedance contrast between water and gas seeping from the seafloor (e.g. Crutchley, Geiger, et al., 2010;Greinert et al., 2010) (Figure S4a). Most of the anomalies, however, are characterised by lower backscattering with a low height to width ratio and disconnection from the seafloor. These appear to result from bio-layers, schools of fish or data acquisition artefacts ( Figure S4b).

| Seismic stratigraphy
We detect 10 unconformities on the seismic reflection profiles based on stratal terminations such as onlap, downlap, toplap and truncation ( Figure 5; Table 3). Unconformities occur as relatively uniform, high-amplitude reflections. They correspond to sequence boundaries U12-U19 documented by Lu and Fulthorpe (2004), as well as two newly identified local unconformities (U12.1, U13.1). The seismic unit above U19 is thickest along the outer shelf and upper slope, and it pinches out at a depth of ca. 500 m on the continental slope.

| Acoustic anomalies
We observe a range of acoustic anomalies in the sub-bottom profiles and MCS reflection data that we summarise in Figure 6 and map spatially in Figure 7. A description of these different acoustic features is given below: 1. Acoustic blanking: Acoustic blanking is a term commonly used in the literature to describe a pronounced reduction in reflectivity due to attenuation of the acoustic signal in the presence of fluid (e.g. Plets et al., 2007;Schroot & Schüttenhelm, 2003). In Figure  6a 7b). These structures can reach depths of 110-800 m and widths of 300 m. The majority of these structures are rooted below U16, whereas their tops coincide with U19. 2. Reflection 'pull-up': Reflection 'pull-up' is a term commonly used to explain the anomalous up-bending of seismic reflections caused by a localised high seismic velocity anomaly (e.g. as seen beneath seabed carbonate mounds in Hovland et al. (1994)). We observe examples of such reflection pull up in MCS data in regions that are up to 500 m wide (Figures 6e and 7c). The most prominent examples of these reflection pull-up zones are capped by high amplitude, positive-polarity reflections (Figure 6e). 3. Reflection 'push-down': In the opposite sense to reflection pull-up features described in (ii), reflection pushdown zones are caused by anomalously low seismic velocities in the overlying sediments (e.g. Arts et al., 2004) (Figure 6d,g). In Figure 6g, we show an example of a push down effect, directly beneath a part of the seafloor reflection that has a distinct reduction in amplitude. The push-down features we observe are clustered along the shelf break, in deeper waters than the pull-up features (Figure 7c). 4. High amplitude reflections: Primarily observed on the MCS data, anomalously high-amplitude reflections extend laterally up to several hundreds of metres. We interpret the polarity of high-amplitude reflections based on the dominant part of the waveform (a trough or a peak) and how that compares to the waveform of the seafloor reflection (see enlarged insets in Figure 6e,f). High amplitude reflections with positive polarity are more common beneath the inner shelf and tend to be associated with the zones of F I G U R E 5 (a) Uninterpreted, and (b) interpreted MCS reflection profile EW00-01-66.  (Figures 6e and 7c,d). Reflections with negative polarity are clustered around the outer shelf to shelf break area, especially in the vicinity of the deep structures with acoustic blanking described in (i) (Figures 6f and 7c,d). Their lateral extent varies from a few tens of metres to >1 km. In the sub-bottom profiler data, anomalously high amplitude reflections are, in places, underlain by broad zones of acoustic blanking, often several kilometres wide (Figures 6h  and 7a,d).

| Sediment geochemistry
The sampled sediments are primarily greenish-greyish muds, infrequently interbedded with thin layers of very fine sand. The ratios of Mn/Al, Ca/Ti and Sr/Ca are presented in Figure S5. Mn/Al ratios are highest in core 29 (at 80 cm and 150 cm below seafloor (bsf)), and cores 26 and 31 (at 20 cm bsf). The ratios are lowest in core 28. The highest Ca/Ti ratios are observed in core 28, followed by cores 30 (at 0, 140 and 250 cm bsf) and 8. The lowest ratios are documented in core 29. Sr/Ca ratios are highest in cores 29, 30 and 31 (at 150 cm, 80 cm and 20 cm bsf respectively), and lowest in cores 8 and 28.

| Pore water chemistry
The pore water chemical measurements are shown in Figure 8. Cl − concentrations are relatively similar for all cores, with the highest values recorded in cores 26 and 31. In core 28, Cl − concentrations show a prominent dip at 330 cm bsf, which is also visible in the Na + concentration. SO 4 2− concentrations are relatively constant between 2700 and 2900 g m − ³ for all cores, except for cores 28 and 31, where they decrease with depth at >150-300 cm bsf. The sharpest decline in SO 4 2− concentrations is observed in core 28 (minimum of 1170 g m − ³). The trends for Ca 2+ concentrations are similar to those of SO 4 2− . A plot of δ 2 H versus δ 18 O shows that pore water values are clustered around the meteoric water line, a reference line that extends between the isotopic compositions of groundwater and seawater (Figure 9). Pore waters from the slope cores (28,29,30,31) exhibit the lowest values of both δ 18 O and δ 2 H, whereas the shelf cores (8, 26) have higher values. Measurements of δ 18 O and δ 2 H ranges between −0.3 and 0.4‰, and −3 and 3‰ respectively. Values of δ 13 C DIC tend to increase with depth (except for core 29) ( Figure  8). Cores 28, 30 and 31 display strong negative deviations, reaching values of −23‰. DIC concentrations display a reverse pattern relative to that of δ 13 C DIC , with the highest values recorded in core 28 at 400 cm bsf.

| Groundwater flow
A comparison of the present-day computed hydraulic head contour maps for all four scenarios (described in section 3.2.3) is presented in Figure 10. We extracted the variation of fluid pressure and hydraulic head with depth at four locations moving from the outer shelf to lower slope (positions A to D) in Figures 11 and 12. The two shoreward profiles represent sections with fine sand and silt facies (positions A and B) while the two distal profiles represent sections entirely comprised of lower permeability clays F I G U R E 7 Maps of acoustic anomalies identified on sub-bottom and MCS profiles: (a) shallow (<50 m) acoustic blanking, (b) deep (>50 m) acoustic blanking, (c) vertical regions of reflection 'pull-up' and 'push-down', (d) high amplitude reflections (positive and negative polarities) (positions C and D) (Figure 11). For the 'high permeability' scenario (scenario 1) with sea-level and sediment loading imposed (Figure 10a), hydraulic heads reach about 1500 m on the continental slope at depth in the thickest region of the model domain, where the strata are clay dominated. Formation pressures reach lithostatic levels at shallow depths in the two distal profiles near the sediment-water interface (Figure 11g,j). However, when silt or fine sand facies are present at the simulated sediment-water interface, fluid pressures are nearly hydrostatic (Figure 11a,d). In all profiles, fluid pressures rise to 40%-60% of lithostatic levels (Figure 11a,d,g,j). For the 'low permeability' scenario (scenario 2; Figure 11b,e,h,k), fluid pressures reach lithostatic levels on the continental slope at nearly all depths for the distal profiles (Figure 11h,k). For the two shoreward profiles, hydrostatic heads still occur at shallow depth.
Scenario 3, in which sea-level loading was removed, is nearly indistinguishable from scenario 2 (compare Figure  10b and 10c; Figure 11b, e, h, k and 11c, f, i, l). When sediment loading is removed (scenario 4), present-day simulated offshore heads are nearly hydrostatic ( Figure  12). Simulated heads for this fourth scenario were as low as 107 m bsl within the clay facies due to the imposed changes in sea level. Because the simulated heads in the fourth scenario were nearly hydrostatic, we did not include them in Figure 11.
Overpressure generation is controlled by the ratio of permeability to specific storage. The vertical permeability of the clay-rich facies that we used in our model is relatively low (10 −19.7 m 2 ), but not unrealistic. Typical measured clay permeability ranges from about 10 −17 to 10 −22 m 2 (Neuzil, 1994). If we had used a higher vertical clay permeability (e.g. 10 −18 m 2 ), lithostatic levels would have not been reached at shallow sub-seafloor depths on the continental slope. Our results are not unique. The head anomalies would have been lower if we had run models with a thinner total sedimentary thickness, because the latter yield a shorter response time of <10 5 years.

| Gas hydrate stability
In the Canterbury Margin, gas hydrates forming from 100% methane are expected to be stable in a minimum F I G U R E 8 Profiles of δ 18 O, δ 2 H, δ 13 C DIC , DIC, Na + , Cl − , SO 4 2− and Ca + concentrations in pore water from sediment cores TAN1703-8, 26, 28, 29, 30 and 31. Locations are shown in Figure 2 water depth of 551 m at a bottom water temperature of 6°C. Our estimated upper limit of methane hydrate stability in sediments is shown by the green contour in Figure  3a. CO 2 hydrates could be stable in minimum water depths of 400 m, at a bottom water temperature of 7°C.
Our estimated upper limit of CO 2 hydrate stability in sediments is shown by the blue contour in Figure 3a. There is good alignment between the upper limit of the slope depressions and up-dip limits of both the methane and carbon dioxide hydrate stability zones.

| Evidence for fluid occurrence, flow and escape
The high amplitude reflections with negative polarity close to the shelf break (Figures 6f and 7d), and the broad zones of acoustic blanking in sub-bottom profiler data on the slope (Figure 7a) are evidence for the occurrence of free gas in sub-seafloor sediments (Judd & Hovland, 1992. Vertical zones of acoustic blanking and reflection 'push-down' (Figures 6d,g and 7a-c), which are often referred to as 'chimneys' or 'pipes', are indicative of focused vertical fluid flow (Cartwright et al., 2007;Løseth et al., 2011). The distinct reduction in seafloor reflectivity above the pushdown effect in Figure 6g is likely caused by free gas in the pore space of near-seafloor sediments. Such reductions in seafloor reflectivity at gas vent locations have been observed elsewhere, including on the Hikurangi Margin, where free gas migrating through the seafloor has reduced the impedance contrast with the overlying water column . The seafloor depressions (Figures 3 and 4), on the other hand, are interpreted as pockmarks, based on their shape and occasional spatial association with chimneys. The gas flares ( Figure S4) provide evidence for localised fluid escape. Localised acoustic blanking features (often with funnel-shaped expressions) in the sub-bottom profiles (Figure 6a-c) are similar to features reported by Hill et al. (2004) and Fulthorpe and Austin (2004) in the US Atlantic margin, and by Stott et al. (2019) in the Chatham Rise. Hill et al. (2004) have interpreted the features as sites where the internal stratigraphic layering was physically disturbed by fluid expulsion, whereas Stott et al. (2019) postulated that the features were buried pockmarks formed at glacial terminations. Fulthorpe and Austin (2004) interpreted the features as associated with catastrophic erosion and redeposition following multiple breaching of glacial lake dams at 12-19 ka. There is no evidence of similar events along the eastern coast of the South Island of New Zealand. We, therefore, consider that the interpretations of Hill et al. (2004) and Stott et al. (2019) better account for the formation of acoustic blanking features in the Canterbury Bight. Our inference is primarily based on the features having similar shapes and sizes to the pockmarks occurring at the seafloor, as well as the spatial correlation to pockmarks where the features intersect the seafloor.

F I G U R E 9
Cross plot of δ 18 O and δ 2 H from cores TAN1703-8, 26, 28, 29, 30 and 31. Locations are shown in Figure 2. The meteoric line is calculated according to Stewart and Taylor (1981) Other anomalies in the seismic reflection data are indicative of sedimentological variability, rather than the presence of fluids. We interpret the anomalously strong, positive-polarity reflections (Figures 6e and 7d) as the seismic response to the occurrence of gravelly and sandy layers extending from onshore, or coarse terrigenous infill of buried palaeo-channels (Bostock et al., 2019;Browne & Naish, 2003;Herzer, 1981;Hillman et al., 2017;Micallef et al., 2020). We interpret the pull-up features (Figures  6e and 7c) as local velocity artefacts due to the shallower occurrence of coarse-grained (high seismic velocity) material. Broad zones of acoustic blanking in sub-bottom profiler data on the shelf (Figure 6h), associated with the strong, positive-polarity reflections, are likely the result of pronounced signal attenuation caused by the presence of overlying coarse-grained (high acoustic impedance) sediment (Stevenson et al., 2002).

| Methane
The sulphate-methane transition zone is defined as the sedimentary interval characterised by a mutual depletion of methane and sulphate due to microbial anaerobic oxidation (Whiticar & Faber, 1986). The pore water chemical record of cores 28, 30 and 31, located on the slope, allow the inference of a sulphate-methane transition zone (Figure 8). The decline in pore water sulphate concentrations and δ 13 C values, associated with an increase in DIC concentrations, suggests that the depth of the sulphatemethane transition zone is likely at ca. 10 m. The associated overall decrease in Cl − and Na + and increase in δ 18 O values, combined with the low Mn/Al ratios in cores 28, 30 and 31, are evidence of the anaerobic oxidation of methane (Bohrmann & Torres, 2006;Jaccard et al., 2009;Li et al., 2016). The shelf cores are too short to allow identification of a similar chemical signature, although a similar decrease in δ 13 C values in the longer core could be postulated (Figure 8).
Methane concentrations measured at scientific drill sites on the outer shelf and upper continental slope ( Figure  2a; IODP Expedition 317 Sites U1351 and U1352, at water depths of 122 m and 343 m respectively), reach 10,000 parts per million by volume (ppmv) and 100,000 ppmv, respectively, in the first 100 m bsf (Fulthorpe et al., 2011). At these two sites, the methane/ethane ratio (C1/C2) is unusually low (<1000 for U1351 and ca. 16,000 for U1352) (Fulthorpe et al., 2011), and the sulphate-methane transition zone is at 15-17 m bsf and lies within the unit above U19, which represents the Last Glacial Maximum (Table  3). At both sites U1351 and U1352, heavier hydrocarbons (C2-C3-C4-C5) were reported in increasing amounts with depth (Fulthorpe et al., 2011). Site U1352 also shows the presence of CO 2 in concentrations of <3000 ppmv (Fulthorpe et al., 2011). At ODP Leg 181 Site 1119 ( Figure  2a; at 393 m water depth), methane values are high again, reaching a concentration of 35000 ppmv in the first 100 m (Carter et al., 1999). The sulphate-methane transition zone is found at 20 m bsf. Values of C1/C2 ratio at Site 1119 are high in all samples. Here, neither heavier hydrocarbons nor CO 2 were detected.
The above geochemical observations indicate that methane is the main fluid forming flow and escape features across the Canterbury Margin. Core data suggest that a biogenic source of methane is more likely at Site 1119. Biogenic gas production is expected when sedimentation rate vs. geothermal gradient yields heating rate values of 7-18°C/Ma, and when organic matter content is at least >0.2% (Clayton, 1992;Torelli et al., 2021;Wellsbury et al., 2002). Across the Canterbury shelf, sedimentation rates F I G U R E 1 0 Present-day computed heads (m) for four scenarios: (a) scenario 1 (high k, sediment and sea-level loading), (b) scenario 2 (low k, sediment and sea-level loading), (c) scenario 3 (low k, sediment loading), (d) scenario 4 (low k, sea-level loading). The four vertical lines denote the positions of hydraulic head and pressure profiles extracted from these simulations estimated from the upper units range between 150 and 250 m/Ma (Lu & Fulthorpe, 2004;McHugh et al., 2018). Based on a geothermal gradient of 36-46°C/km, and an organic matter content of 0.6%-1% for both sites U1351 and U1352 (Fulthorpe et al., 2011), we can deduce that conditions were conducive to biogenic gas production across the shelf. On the other hand, the vertical variability of heavier hydrocarbons at sites U1351 and U1352, the low C1/C2 ratios, the deep and the presence of CO 2 at site U1352 lead us to also infer a deeper source of methane, likely with an early thermogenic origin (Fulthorpe et al., 2011;Milkov & Etiope, 2018). Chimneys/pipes provide the main conduits for methane from the deeper reservoirs (Figures 6d and  7b). Previous studies (based on seismic facies interpretation  Figure 10 2015). Without a chemical analysis of the seeping fluids, it is difficult to determine whether biogenic and/or thermogenic methane is responsible for pockmark formation.

| Groundwater
Borehole records from across the Canterbury Margin generally document pore waters with salinities equivalent to seawater down to ca. 300-500 m bsf (Carter et al., 1999;Fulthorpe et al., 2011). However, one borehole (IODP Expedition 317 Site U1353) shows fresher pore water salinities (24 practical salinity units) down to depths of 180 m bsf (Expedition 317 Scientists, 2011). Recently acquired seismic reflection and electromagnetic data point to the occurrence of an extensive offshore freshened groundwater system across the Canterbury Bight, with a maximum thickness of at least 250 m and a top that reaches a minimum depth of 50 m bsf (Micallef et al., 2020). The offshore freshened groundwater system extends at least 60 km perpendicularly from the coastline. However, the offshore extent may be higher because pore water records from Site 1119 still show some freshening immediately below the seafloor at a distance of ca. 90 km from the coastline (Carter et al., 1999). Based on numerical modelling, Micallef et al. (2020) inferred that the groundwater system was predominantly recharged during sea-level lowstands. At these times, topographically-driven flow took place across much of the continental shelf, and groundwater velocities were about an order of magnitude higher than during sea-level highstands. The δ 18 O and δ 2 H isotopic depletion of pore waters on the slope, relative to those on the shelf (Figure 8), suggests that a mixture of saline and freshened groundwaters are being expelled across the slope.

Sediment loading
Our numerical models suggest that sediment loading can generate pore pressures that reach lithostatic levels across the continental slope (Figures 10 and 11). Such pore pressures are sufficient to entrain and disperse seafloor sediments and form pockmarks in both the high and low permeability scenarios. At shallower water depths, pore pressures only reached lithostatic levels in the outer shelf (position B) for the low permeability scenarios where clay sediments are present (Figures 10 and 11). The latter is likely a result of low sedimentation rates due to low accommodation space. In addition, during the last glacial cycle, sedimentation rates appear to have been higher on the outer shelf and upper slope (140-500 m water depth), and lower elsewhere (unit above U19 in Figure 5a,b). Interestingly, the seafloor pockmarks are located where the unit above U19 is thinnest or absent (Figure 5a,b).
The above leads us to infer that: (i) overpressure generated by sediment loading more likely accounts for the formation of pockmarks on the continental slope than on the shelf; and (ii) the weight of overburden above the deeply buried and overpressured interval in the outer shelf (position B in Figures 10 and 11) is thwarting fluid expulsion at the seafloor in this location.

Sea-level loading and offshore freshened groundwater flow
The hydraulic heads associated with sea-level loading and the emplacement and flow of offshore freshened groundwater, particularly at sea-level lowstands, are never high F I G U R E 1 2 Computed present-day head profiles for the four scenarios (1-4), at the four selected positions A, B, C and D. Positions are shown in Figure 10 enough in our models to independently reach lithostatic levels across either the shelf or slope. However, they may play a minor role by supplementing fluid pressures associated with sediment loading and methane generation.

Gas hydrate dissociation
Across the continental slope, we observe that the shallowest slope pockmarks roughly coincide with the upper limit of methane and carbon dioxide hydrate stability (Figure 3). We rule out the role of gas hydrate dissociation, however, because we do not observe BSRs in any of the available seismic data, even in the deeper sections of the margin. If gas hydrate dissociation had occurred as a result of sea-level lowering and/or an increase in bottom water temperatures, we would still expect to see BSRs (or other seismic manifestations of a gas hydrate system) in deeper waters (i.e. >700 m bsl). In addition, the porewater geochemical data from Site 1119 show that the variations of Na + and Cl − with depth in the upper 100 m are dissimilar and they do not replicate the variation of salinity with depth (Carter et al., 1999). This suggests that methane hydrate dissociation did not play a role in pore water freshening at Site 1119 in the last 250 ka.

Other processes
The processes listed above do not generate the overpressures required to form pockmarks across the shelf. Other processes associated with methane generation and flow (e.g. reservoir overcharging, buoyancy) (Hedberg, 1974;Law & Dickinson, 1985;Timko & Fertl, 1971;Vargas-Silva et al., 2019;Zhao et al., 2018) as well as a reduction in hydrostatic pressure due to sea-level lowering (Riboulot et al., 2014), could have played a role. We do not have enough information at the Canterbury Margin to numerically test these hypotheses.

| Sediment type
The distribution of the pockmarks appears to be influenced by sediment distribution across the seafloor. Pockmarks reach densities of 100 per km 2 in the muddy areas of the outer shelf and lower slope, whereas they are absent in the sandy areas of the inner shelf (Bostock et al., 2019;Browne & Naish, 2003;Herzer, 1981) (Figures 1b and 3a). Pockmarks are less likely to form in sands because these are difficult to disperse and allow faster flow of fluids (e.g. Chand et al., 2012).

| Bottom currents
The sedimentation rate and pockmark morphology also appear to be regulated by bottom currents. On the slope, the west-south-west to east-north-east asymmetry of the pockmarks at <500 m bsl (Figures 3a and 4b), and possibly the thin/absent U19-seafloor sediment package (Figure 5b), are a result of the influence of a north-east flowing bottom current in the slope area. The Southland Current, STF and Antarctic Intermediate water all flow in this direction at velocities of ca. 20 cm s −1 (Chiswell, 2003;Hadfield et al., 2007), which is sufficient to transport or erode silt to clay sediments (Hjulstrom, 1939). The currents were likely stronger during glacial periods due to strengthening of temperature gradients and wind regimes (Davy et al., 2010). Most shelf pockmarks have a dominant north-east to south-west orientation and an inverse asymmetric profile relative to those on the slope (Figures 3a and 4a). This geometry suggests interaction with a southward-flowing bottom current, but no such southward flowing water mass has been identified to date. However, water temperature differences between the southern and northern Canterbury Bight in the 100-200 m depth range suggest the occurrence of turbulence in the Southland Current, which could create eddy structures inshore that may be responsible for the asymmetry in pockmark morphology (Beentjes et al., 2002).

| Timing of fluid expulsion
Evidence for present day seafloor fluid seepage is scant and includes <10 flares over the entire study area. Only ca. 5% of the shelf pockmarks have the circular planforms and higher backscatter than the surrounding seafloor, which are indicative of recent formation due to insufficient time to develop extensive elongation by bottom currents or burial by sediments. The occurrence of the sulphate-methane transition zone at depths of >10 m confirms that there is no methane seepage at the seafloor at the sampled sites at present (e.g. Bhatnagar et al., 2008;Borowski et al., 1996;Coffin et al., 2008;Mogollon et al., 2012); this is corroborated by the fact that the top of most chimneys/pipes coincides with U19. Pockmarks formed by methane seepage are thus likely to have developed in the past. Pockmarks associated with overpressured saline and freshened groundwaters across the slope, on the other hand, may be relatively younger, because pore pressures due to sediment loading reach lithostatic levels at present (Figures 10 and 11).

| CONCLUSIONS
We have investigated pockmarks across the Canterbury Margin-using multibeam echosounder and seismic reflection data, sediment cores and pore water samples and numerical models of groundwater and gas hydrates-to identify the drivers and controls of their formation. Our main conclusions are the following: 1. There are >6800 surface pockmarks with densities of up to 100 per km 2 , and an undefined number of buried pockmarks, in the middle to outer shelf and lower continental slope of the Canterbury Margin. 2. Geophysical evidence for the occurrence of free gas in the sub-seafloor includes high amplitude reflections with negative polarity close to the shelf break and broad zones of acoustic blanking in sub-bottom profiler data on the slope. Shallow to deep focused vertical fluid flow in the sub-seafloor is indicated by vertical zones of acoustic blanking and reflection 'push-down' (chimneys/pipes). 3. Methane-biogenic and/or thermogenic in origin-is the main fluid forming flow and escape features across the Canterbury Margin, although a mixture of saline and freshened groundwaters may also be seeping across the slope. 4. The drivers of fluid flow and seepage include overpressure across the continental slope due to sediment loading and thin sediment overburden above the overpressured interval in the outer shelf. Processes associated with methane generation and flow, and a reduction in hydrostatic pressure due to sea-level lowering, may also account for fluid flow and seepage features, particularly across the shelf. 5. Pockmark formation is restricted to regions of the outer shelf and slope where there is a higher concentration of muds at the seafloor. 6. After formation, pockmarks across the slope have been elongated by bottom currents (e.g. Southland Current, STF and the Antarctic Intermediate Water flow). Elongation of shelf pockmarks may be related to eddy structures of the Southland Current. 7. Whereas pockmarks formed by methane seepage developed in the past, pockmarks arising from the seepage of overpressured saline and freshened groundwaters across the slope could have formed at present.
In comparison to previous studies of the Chatham Rise and Otago Margin, the interdisciplinary approach employed in our study has demonstrated that, in the Canterbury Margin, pockmark occurrence, characteristics and formation can be influenced by a combination of multiple factors that include different fluid types, overpressures, seafloor sediment type, stratigraphy and bottom currents.