Origin of High Density Seabed Pockmark Fields and Their Use in Inferring Bottom Currents

: Some of the highest density pockmark ﬁelds in the world have been observed on the northwest Australian continental shelf (>700/km 2 ) where they occur in muddy, organic-rich sediment around carbonate banks and paleochannels. Here we developed a semi-automated method to map and quantify the form and density of these pockmark ﬁelds (~220,000 pockmarks) and characterise their geochemical, sedimentological and biological properties to provide insight into their formative processes. These data indicate that pockmarks formed due to the release of gas derived from the breakdown of near-surface organic material, with gas accumulation aided by the sealing properties of the sediments. Sources of organic matter include adjacent carbonate banks and buried paleochannels. Polychaetes biodiversity appears to be affected negatively by the conditions surrounding dense pockmark ﬁelds since higher biodiversity is associated with low density ﬁelds. While regional bi-directionality of pockmark scours corresponds to modelled tidal ﬂow, localised scattering around banks suggests turbulence. This multi-scale information therefore suggests that pockmark scours can act as proxy for bottom currents, which could help to inform modelling of benthic biodiversity patterns. area. Similar results were observed for the scoured pockmark.


Introduction
Pockmarks are seafloor depressions that form in soft sediment by fluid-flow processes. They have been reported since the 1970s [1] and are present worldwide [2]. Pockmarks are most commonly found in areas of thick Holocene sediment accumulation along continental shelves, in tropical to glacial settings. The expression of fluid-flow via pockmarks at the seafloor is important for several reasons. These include the need to identify potential hydrocarbon resources and areas where carbon sequestration is not possible due to the potential leakage of the targeted reservoir [3]; areas of potential geohazards (e.g., gas venting, slope failure, basement faulting) [4,5]; and areas of significant biodiversity associated with seep-related habitats, such as bioherms and reefs [1,6,7].
(1) Identify and map pockmarks on the shelf of the western Timor Sea, Australia using a new semi-automated method. (2) Document the geochemical and sedimentological properties of pockmark sediments, and of associated infaunal communities. (3) Present a conceptual model for pockmark formation that links pockmark form and density to local environmental conditions, including hydrodynamics and sedimentary processes.

Geological and Physiographic Settings
The study area is located in the tropical Oceanic Shoals Australian Marine Park (AMP; 71,740 km 2 ) within the Timor Sea on the continental shelf of northwestern Australia. The area lies within the sedimentary Bonaparte Basin, which is expressed at the seabed by the Van Diemen Rise, Sahul Shelf, Malita Graben, Sahul Syncline, Petrel Sub-basin and the Londonderry Rise ( Figure 1). The Sahul Shelf, Van Diemen Rise and Londonderry Rise form carbonate-dominated platforms characterised by a complex suite of submerged flat-topped carbonate banks and terraces, separated by valleys, channels and plains [21][22][23][24][25]. Water depths range from 10 m on the shallowest banks to over 200 m in shelf-incising valleys. The central basin is about 270 km long with water depths of 90 to 140 m. The seabed of the greater Bonaparte Basin, including the Oceanic Shoals AMP region represents a drowned landscape consisting of coastal lowlands, estuarine valleys and nearshore environments. These are the product of periodic and repeated formation of a semi-enclosed sea during the Quaternary and most recently inundated during the Late Pleistocene to Holocene postglacial rise in sea level [26][27][28] Australian Marine Park (brown outline) in the Timor Sea, Northern Australia. Geological setting is also included with the encompassing Bonaparte Basin as a bold dashed line and its subdivisions as continuous line. Post-survey reports that are associated with the survey outlined in this figure are as followed: Oceanic Shoals [21], Petrel Sub-basin [29], LaPerouse-Durville [29], Darwin Shelf [24,25], cores [27]. Figure modified from [21] 1. 1

.2. Environmental and Oceanographic Settings
The Oceanic Shoals study area is within a tropical region, with a pronounced southern hemisphere summer monsoonal season that commonly includes cyclones. The Indonesian Through-Flow, passing between Timor and northern Australia delivers relatively warm (~26 • C to 31 • C) low salinity water from the north. Across the shelf, macrotidal conditions drive strong tidal flows and wind-driven waves and seasonal cyclones come from the east. These conditions result in an anticlockwise ocean circulation across the shelf strongly influencing sediment transport [30]. The re-suspended seabed sediments coupled with sediment discharged from coastal river systems, contribute to a turbid shelf environment that is likely to suppress photosynthetic activity at the seabed. This therefore may be an important determinant of the composition of benthic assemblages across the region [31].
The study area encompasses three areas located along the north and west portion of the central basin ( Figure 1). Of these, Area 1 is located within the Petrel Sub-basin, Area 2 straddles the southern portion of the cross-shelf valley within the Sahul Syncline (Figure 1), and Area 3 is located across the Sahul Syncline to Malita Graben transition. The extent of this general area spans a cross-shelf transect of approximately 110 km ( Figure 1). Each area was positioned to encompass a variety of seabed geomorphic features and water depths, as identified in the national geomorphic features database [21,23]. On most of the Australian continental shelf, pockmarks are seldom found because of the dominance of coarse sediments (sand and gravel) in most regions. There is also a limited input of terrigenous organic-rich sediments, and a moderate to high percentage of carbonate sediments on the shelves [27,30]. However, high-density fields of pockmarks are present within the broad continental shelf of northern Australia, including the extensive Bonaparte Basin where seabed sediments do incorporate terrigenous fine-grained material. During the Quaternary, the present seabed of the Bonaparte Basin has been repeatedly exposed either in part, or completely, from eustatic sea level changes [27,28,32]. During low sea-levels, banks, shoals and terraces became exposed, resulting in Joseph Bonaparte Gulf forming a semi-enclosed basin [33]. Within the Bonaparte Basin, the Petrel Sub-basin, Malita Graben, Sahul Syncline and adjacent structural elements have acted as a depocentre for marine and riverine sediment since the Paleozoic [26,27,[33][34][35][36]. This resulted in hydrocarbon generation occurring at depths in older sedimentary strata. Recently, seabed sediment with a particularly high organic matter content derived from the preservation of mangroves and other plant matter has been reported within paleo-river channels mapped on the Petrel Sub-basin mid-shelf and presently covered by a thin sediment veneer [33]. Additionally, the reworking of buried terrestrial sediments upwards to the seabed via pockmark formation has been identified within the Petrel Sub-basin using rare earth element chemistry [37]. Previous studies suggested that banks of the greater Oceanic Shoal region could result from one of these two processes: (1) hydrocarbon release from depth and its movement upward in faults providing a food source for the colonising reef-builders, or (2) thinner sediment cover over organics resulting in faster release of the product of organic breakdown [29].

Survey Summary
The data presented herein were collected during the Oceanic Shoals environmental marine survey (SOL5650/GA0339) undertaken in September 2012 [21]. High-resolution multibeam echosounder (MBES) bathymetry and backscatter data were collected over three areas covering a total of about 425 km 2 of seabed ( Figure 1). MBES bathymetric data was processed using Caris Hips and SIPS v.7.1. Seabed backscatter data was processed with CMST-GA MB Process toolbox software co-developed by the Centre for Marine Science and Technology at Curtin University of Technology and Geoscience Australia [38]. Both data types were gridded at 1 m horizontal resolution for analysis and interpretation purposes.
Shallow sub-bottom profiles (SBP) were collected over 140 line km using a sparker source and processed with a C-view processing unit using 1 sec firing rate and 500 ms record length, then filtered using band-pass filter and converted to SGY file. Seabed samples were collected in duplicates at 61 stations using a Smith-McIntyre grab or box corer ( Figure 2). Each sample was analysed for sedimentological characteristics (sediment texture as Mud, Sand and Gravel %, grain-size distribution, CaCO 3 content) and geochemical composition (incl. porosity, CO 2 production rates (TCO 2 ), and Total Organic Content (TOC)) as described in [21,39] method, and polychaetes biodiversity (57 stations only-identified to species level according to [17]).

Geomorphic Analysis and Automated Pockmarks Characterisation
Geomorphic feature mapping was completed manually using MBES data and derivative datasets such as slope and contour to identify discrete features ( Figure 3) [21]. Feature definitions were based on [23]. For this study, bank morphometric profiles were calculated using ArcGIS 3D (v. 10.1, Esri, Redlands, CA, USA) profiling across the narrowest axis of the bank (Figure 3b).  Due to their large number, manual mapping of individual pockmarks was deemed unsuitable. Therefore, an automated method (GA1) was developed to identify both pockmarks and associated scour marks. The overall approach was similar to other published automated and semi-automated methods [10,18], in which a multi-stage process was undertaken. The main stages included (1) identifying pockmarks and differentiating scoured from non-scoured pockmarks, (2) determining scour orientation, (3) calculating field density and (4) assessing confidence in the results ( Figure 4). Each step of the method used an existing suite of ESRI spatial analytical tools, such as the hydrological and surface toolboxes. Details of the method can be found in the supplementary material. Density classes were used as a common unit for comparative analyses with other survey datasets. Due to their large number, manual mapping of individual pockmarks was deemed unsuitable. Therefore, an automated method (GA1) was developed to identify both pockmarks and associated scour marks. The overall approach was similar to other published automated and semi-automated methods [10,18], in which a multi-stage process was undertaken. The main stages included (1) identifying pockmarks and differentiating scoured from non-scoured pockmarks, (2) determining scour orientation, (3) calculating field density and (4) assessing confidence in the results ( Figure 4). Each step of the method used an existing suite of ESRI spatial analytical tools, such as the hydrological and surface toolboxes. Details of the method can be found in the supplementary material. Density classes were used as a common unit for comparative analyses with other survey datasets. The pockmarks identified using the GA1 method were compared with results from two other methods using a small subset area of Area 3. These methods were a combination of the British Geological Survey (BGS) Seabed Mapping toolbox [18,20] and the Geomorphometry and Gradient Metrics (GGM) toolbox [40], and another tool tested by Geoscience Australia (GA2). The subset area included a carbonate bank and the associated clustered pockmarks surrounding its base, as well as other scoured and non-scoured pockmarks distributed in the plain area.
For the BGS combined approach, the landform tool in the GGM toolbox was used to create the surface curvature landform index (Bolstad's variant) using a square neighborhood analysis window of 12 by 12 cells. The derived raster was used as input for the Feature delineation [Derived] tool within the BGS Seabed Mapping Toolbox [20]. The mapped features corresponded to the areas with a surface curvature landform index lower than −0.04, that are wider than 5 m with a minimum size The pockmarks identified using the GA1 method were compared with results from two other methods using a small subset area of Area 3. These methods were a combination of the British Geological Survey (BGS) Seabed Mapping toolbox [18,20] and the Geomorphometry and Gradient Metrics (GGM) toolbox [40], and another tool tested by Geoscience Australia (GA2). The subset area included a carbonate bank and the associated clustered pockmarks surrounding its base, as well as other scoured and non-scoured pockmarks distributed in the plain area.
For the BGS combined approach, the landform tool in the GGM toolbox was used to create the surface curvature landform index (Bolstad's variant) using a square neighborhood analysis window of 12 by 12 cells. The derived raster was used as input for the Feature delineation [Derived] tool within the BGS Seabed Mapping Toolbox [20]. The mapped features corresponded to the areas with a surface curvature landform index lower than −0.04, that are wider than 5 m with a minimum size ratio higher than 0.2 and that include cells with values lower than −0.05. The final outline was achieved by adding a 2 m buffer.
The GA2 method shared similarities with BGS approach, in which Hydrology toolset in ArcGIS TM was used to identified sinks and pits from the bathymetry data. However, GA2 used a fuzzy features tool in landserf (based on topographic terrain analysis) to calculate multiscale pitness [41] and the openness tool (based on the viewshed principle) to calculate positive openness [42]. Both calculations were based on a search radius of 21 cells. Where pitness was greater than the mean plus one standard deviation, pits were identified. From these, only the ones associated with a positive openness values of less than the mean minus one standard deviation were selected and converted into polygons (simplify boundary function ON). Selection was further constrained by only selecting polygons that intersected either sinks, or pits or both. Finally, polygons with area less than 10 m 2 were deleted and resultant polygons were regarded as pockmarks. The key issue related to these semi-automated methods is the choices of value thresholds, which is a somewhat subjective process.

Infaunal Data Analyses on Grab Samples
Multivariate analyses were performed only on samples that returned polychaetes identified to species level, while univariate analyses of species richness and abundance included samples from which no polychaetes were identified. All samples from banks were excluded from statistical analyses because these did not have an associated value for pockmark density. Assemblage data and species richness data were square-root transformed to reduce the influence of dominant groups [43] and meet ANOVA assumptions, respectively. A PERMANOVA using unrestricted permutation of raw data and type III SS [44] was undertaken to determine if polychaete assemblages varied among pockmark density classes, while a 1-way ANOVA was undertaken to determine if species richness varied among pockmark density classes.

Hydrodynamic Model
To understand the processes contributing to the presence of pockmarks and associated scouring, a coastal shelf model was constructed to evaluate bed shear stress at the three different survey areas ( Figure 1). A finite element model utilising an unstructured grid was constructed. The bathymetry used in the model was Australian Bathymetry and Topography Model [45] gridded at 250 m and the multibeam dataset used here gridded at 10 m resolution.
The boundary of the model was set beyond the shelf break in order for the currents to stabilise by the time they reach the study area. The northern boundary conditions were set as open water boundaries between islands of the Indonesian archipelago and Papua New Guinea. Both the eastern and western boundaries were set as open water. Two types of ocean boundary conditions were trialed and both included the semi-diurnal tides. (1) Tidal elevations derived from the TOPEX/Poseidon satellite altimetry global tidal model driver (TMD) toolbox version 2.05 with all constituents implemented to create hourly tidal elevation files, and (2) tidal currents derived from the same toolbox using all harmonic constituents. These currents were multiplied by the width and depth of the boundary to calculate the hourly flow rate. Wind stresses from the Bureau of Meteorology gridded datasets were also applied in both of the seasonally dominant directions of the North West and South East. The wind forcing, while changing the surface currents, did not change the currents near the sea bed and the bed shear stress calculations. The model was refined near the study area and, with salinity kept constant, run for a three-month period to represent mean tidal conditions. From these model runs, the bed shear stress time series and exceedance of erosion was also calculated.

Geomorphic Analysis
High-resolution seafloor mapping revealed multiple carbonate banks rising tens of meters above otherwise flat pockmarked plains [21].

Banks
Relatively steep-sided (slopes up to 35 degrees) and flat-topped banks are located in water depths of 45 m to123 m ( Figure 2). Along their flanks, banks have terraces of various widths, which are regarded as possible indicators of lower sea levels. Banks and their associated features are characterised by the highest backscatter intensity values measured from the multibeam sonar data and coarsest sediment composition sampled [21]. Sub-bottom profiles reveal a mostly acoustically transparent seismic facies across the banks (Figure 2). The dimensions (length and height) of the banks reduces northward towards the shelf edge ( Figure 3).

Plains
Plains dominate the seascape, occur between 94m and 146 m water depth and deepen northward. They record gradients averaging 2.5 degrees and are characterised by low backscatter intensity values, representative of the fine-grained sediment observed in the samples [21]. Their sub-surface is mostly characterised by a thick sequence (>150 ms TWT) of acoustically transparent seismic facies closest to the seafloor, underlain by a semi-stratified facies, which in places is replaced by cut and fill facies (Figure 2d). The thickest section (>350 ms TWT) is recorded in the deepest region of Area 2, which corresponds to the start of a deep cross-shelf valley (Figures 1 and 2).

Pockmarks: Density, Scour Direction Results
Semi-automated mapping identified 220,508 pockmarks throughout the plain areas, including at the foot of the banks. Some bank tops also hosted pockmarks, but these were not included in this analysis. Pockmarks displayed wide ranging geometries, varying from circular to elongate. Their dimensions varied from 10 to 30 m in diameter and 0.5 to 2 m in depth. In most cases, the centers of the pockmarks recorded higher backscatter intensity values than their surroundings [21].
Elongated pockmarks represented mainly pockmarks accompanied by a semi-linear depression shallowing in a particular direction (Figure 3b). The depression, commonly but not always, formed a v-shape on the seabed, with the open end of the 'v' on the downstream side of the pockmark. These depressions extended for up to~200 m, but on average less than 80 m, and are interpreted as a scour mark ( Figure 5). For the purpose of this study, pockmarks are divided in two main types: with and without scour. Pockmark types are intermixed throughout the areas they are found within and pockmarks without scour are at most 4.5 times more common (see supplementary tables). Scoured pockmarks show strong bi-directionality, which is associated with an overall east-west quadrant ratio of~1 ( Figure 5). The scour directions in Areas 1 and 2 are slightly more biased towards the eastern quadrant, whereas the directions recorded in Area 3 dominate the western quadrant. More specifically, the bi-directionality is strongest in Area 3, with over 75% of the scour directions divided relatively equally between W-NW and E-SE quadrants, while Area 1 is the most scattered (44% within the W-NW and E-SE quadrant; Figure 5 and supplementary tables).
Pockmarks form fields of dense clusters, with densities greater than 700/km 2 for non-scoured pockmarks ( Figure 6). High-density areas are mostly concentrated around the foot of the banks and represent over 17% of the total surveyed area (supplementary Figure S2). The most common density for both types of pockmarks is between 25 and 200/km 2 . Area 1 contains the highest coverage of banks and low density pockmark fields (<25/km 2 ; Figure 6 and supplementary Figure S2). Area 2 contains a large field of high density pockmarks within the cross-shelf valley (Figures 2 and 6) and the largest area of intermediate-density (200-700 pockmark/km 2 ) fields of non-scoured pockmarks ( Figure 6 and supplementary material). Many of the scoured pockmarks in the central plain of Area 2, which extends across the southern end of the Sahul Valley (Figure 1), are oriented so that a flow from north to south is indicated for the Sahul Valley. This is consistent with the general anticlockwise circulation suggested for the basin [15]. In Area 3, dense non-scoured pockmarks fields surround the banks and similarly to Area 2, many of the pockmark scours are oriented such that a flow from northwest to southeast is implicated, also consistent with the anticlockwise circulation [30]. 2, which extends across the southern end of the Sahul Valley (Figure 1), are oriented so that a flow from north to south is indicated for the Sahul Valley. This is consistent with the general anticlockwise circulation suggested for the basin [15]. In Area 3, dense non-scoured pockmarks fields surround the banks and similarly to Area 2, many of the pockmark scours are oriented such that a flow from northwest to southeast is implicated, also consistent with the anticlockwise circulation [30].

Comparison Analysis and Confidence Assessment
Comparison analysis between the results of GA1, GA2 and BGS methods shows that all three methods identified most pockmarks present, within 10% of each other (Figure 7). However, discrepancies exist between their pockmark total count and coverage area (Figure 7i). For example, BGS method identifies 60% more pockmarks than the GA1 method (Figure 7i), but their total coverage areas are similar. These discrepancies are attributed to GA methods often combining closely spaced pockmarks or identifying pockmark scours. For GA1 method, this resulted respectively in reducing the overall pockmark count while increasing the total coverage area (Figure 7i).  Panels (a-d) show the picks using the various methods for the entire test area, while panels (e-h) show an inset of the test area in order to present more details. (i) Histogram comparing the total pockmark count (grey bars) and coverage area (coloured bars) between the three automated methods used to identify pockmarks.

Sediment Texture Analyses
Sediment textures provide some insight into pockmark development. Seabed sediments in the mapped areas have mixed textures, including carbonate-dominated sand (90-95% CaCO 3 ) and gravel (~95% CaCO 3 ), and calcareous mud (CaCO 3 values around 25-35%). Some of these sediment samples have the highest mud concentrations measured from the Australian continental shelf [40]. Most of this mud is composed of silt-sized particles, with clay content up to~25% (Figure 8). The clay is expected to be about 50% smectite, based on data from seabed sediments in the Eastern Joseph Bonaparte Gulf [46]. Low-density pockmark fields have the lowest mud content and highest sand content (Figure 9). This relationship reverses for higher density fields but mud content is greatest in the plain areas where pockmark densities are less than 200/km 2 . Gravelly mud and muddy gravel generally occur where the highest pockmark field densities are found (i.e., foot of the banks), with the mud to sand ratio less pronounced than for the plains (Figure 9). These textural relationships are also directly correlated with the average seabed backscatter intensities, in which the highest backscatter intensities are associated with the banks and where the sand to mud ratio is greater than 1. The opposite relationship is observed for the areas of lowest backscatter intensities ( Figure 9).

Geochemistry Analyses
Oceanic Shoals sediments have the highest Total Organic Carbon (TOC) concentrations of sediments collected so far on the Australian continental shelf, which average 0.8% (Figure 10a, [39]). Moreover, there is a positive linear relationship between the mud content, surface area (SSA) and Figure 9. Sediment characteristics (box plots) and averaged seabed backscatter strength (coloured symbols) for each non-scoured pockmark field density class per area. Similar results were observed for the scoured pockmark.

Geochemistry Analyses
Oceanic Shoals sediments have the highest Total Organic Carbon (TOC) concentrations of sediments collected so far on the Australian continental shelf, which average 0.8% (Figure 10a, [39]). Moreover, there is a positive linear relationship between the mud content, surface area (SSA) and TOC as common for Australian Shelf sediments [48]. While most Oceanic Shoals sediments have surface area-normalised TOC concentrations that are in the typical range for the continental shelf [39], the sediments from the dense pockmark fields have less TOC per unit surface area than sediments from the low-density pockmark fields (Figure 10b). TOC as common for Australian Shelf sediments [48]. While most Oceanic Shoals sediments have surface area-normalised TOC concentrations that are in the typical range for the continental shelf [39], the sediments from the dense pockmark fields have less TOC per unit surface area than sediments from the low-density pockmark fields (Figure 10b).  Table 1 in [39]). (b) TOC versus specific surface area of sediments highlighting that sediments from the dense pockmark fields have less TOC per square of sediment than sediments from the low-density pockmark fields.
When organic matter breaks down in marine sediments, dissolved inorganic carbon (DIC) is released to porewaters where it may build up in concentration. Porewater DIC concentrations of seabed sediments on Australia's continental margin (as TCO2 pool) increase approximately exponentially with porosity, and the slope of this relationship is about two-fold higher for pockmarked sediment of the Oceanic Shoals area (Figure 11).  Table 1 in [39]). (b) TOC versus specific surface area of sediments highlighting that sediments from the dense pockmark fields have less TOC per square of sediment than sediments from the low-density pockmark fields.
When organic matter breaks down in marine sediments, dissolved inorganic carbon (DIC) is released to porewaters where it may build up in concentration. Porewater DIC concentrations of seabed sediments on Australia's continental margin (as TCO 2 pool) increase approximately exponentially with porosity, and the slope of this relationship is about two-fold higher for pockmarked sediment of the Oceanic Shoals area (Figure 11).
When organic matter breaks down in marine sediments, dissolved inorganic carbon (DIC) is released to porewaters where it may build up in concentration. Porewater DIC concentrations of seabed sediments on Australia's continental margin (as TCO2 pool) increase approximately exponentially with porosity, and the slope of this relationship is about two-fold higher for pockmarked sediment of the Oceanic Shoals area (Figure 11).  Table 1 in [39] for survey details). The yellow triangles pertain to pockmarked sediment in the Oceanic Shoals AMP survey areas.  Table 1 in [39] for survey details). The yellow triangles pertain to pockmarked sediment in the Oceanic Shoals AMP survey areas.

Polychaete Analysis
A total of 38 grab and box core samples were considered for assemblage data. There were three outliers (one for each area) in assemblage data that were identified by their isolation from the main cluster. Statistical analyses were thus performed with and without these (see supplement Figure S5 for details on the distribution). There were no significant differences in polychaete assemblages for pockmark density without scours (including outliers: Pseudo F = 1.2583, p = 0.051; excluding outliers: Pseudo F = 1.1909, p = 0.107). In contrast, there were significant differences in assemblages among pockmark density classes with scours (including outliers: Pseudo F = 1.4036, p = 0.023; excluding outliers: Pseudo F = 1.14548, p = 0.018). Pairwise comparisons revealed significant assemblage (p = 0.005) differences between only the lowest densities classes of scoured pockmarks (<25/km 2 and 25 to 200/km 2 ).
A total of 51 samples were considered for species richness data. Among pockmark density classes without scours, species richness was significantly different (df = 3, F = 8.7349, p = 0.0001). Tukeys HSD pairwise tests showed significantly more species in the lowest density class (<25 pockmarks/km 2 ) than in the other three higher density classes (Figure 12a). Similarly, species richness was significantly different among pockmark density classes with scours (df = 2, F = 8.6768, p = 0.0006), with more species collected in the lowest density class than the higher density classes (Figure 12b).

Polychaete Analysis
A total of 38 grab and box core samples were considered for assemblage data. There were three outliers (one for each area) in assemblage data that were identified by their isolation from the main cluster. Statistical analyses were thus performed with and without these (see supplement Figure S5 for details on the distribution). There were no significant differences in polychaete assemblages for pockmark density without scours (including outliers: Pseudo F = 1.2583, p = 0.051; excluding outliers: Pseudo F = 1.1909, p = 0.107). In contrast, there were significant differences in assemblages among pockmark density classes with scours (including outliers: Pseudo F = 1.4036, p = 0.023; excluding outliers: Pseudo F = 1.14548, p = 0.018). Pairwise comparisons revealed significant assemblage (p = 0.005) differences between only the lowest densities classes of scoured pockmarks (<25/km 2 and 25 to 200/km 2 ).
A total of 51 samples were considered for species richness data. Among pockmark density classes without scours, species richness was significantly different (df = 3, F = 8.7349, p = 0.0001). Tukeys HSD pairwise tests showed significantly more species in the lowest density class (<25 pockmarks/km 2 ) than in the other three higher density classes (Figure 12a). Similarly, species richness was significantly different among pockmark density classes with scours (df = 2, F = 8.6768, p = 0.0006), with more species collected in the lowest density class than the higher density classes (Figure 12b).

Hydrodynamic Model Analysis
The performance of the model developed here was considered good as the model replicated the growth of the tidal amplitude from the edge of the shelf to the shoreline. The model results were compared to the Timor Sea mooring array, part of the Integrated Marine Observation System (IMOS), that runs to the shelf edge and records a tidal amplitude that grows from plus 3 m at the shelf edge to 6 m on the inner shelf. Results were also correlated with other mooring sites and tide gauges, including Darwin Harbour, Melville and Milner Bay, and inshore the Gulf of Carpentaria.
There is no data available on the bed shear stress erosion threshold for this area above which sediment transport would occur. However, a value of 0.1 N/m 2 is considered an appropriate value for reasonably cohesive sediments, such as the ones found in this study area [49] Over the three months period modelled for the study area, bed shear stress exceeded 0.1 N/m 2 for over 33% of the time and 0.2 N/m 2 for over 8% of the period (Figure 13). Area 1 is exposed to the longest period of bed shear stress exceeding the threshold (49%), while Area 2 the least (12%). Over periods of many years, even the periods of high bed shear stress would be sufficient to provide scour energy on the seabed. Analyses of the modelled current velocities (directions) for all areas suggest a strong bi-directionality with a slight dominance of the current in the NW direction ( Figure 14).

Confidence Assessment of the Semi-Automated Method
Method GA1 used to identify pockmarks and accompanying scour mark was compared to the BGS combined delineation approach and a second unpublished method developed by Geoscience Australia (GA2). The comparative analysis shows that while all pockmarks in the area were

Confidence Assessment of the Semi-Automated Method
Method GA1 used to identify pockmarks and accompanying scour mark was compared to the BGS combined delineation approach and a second unpublished method developed by Geoscience Australia (GA2). The comparative analysis shows that while all pockmarks in the area were identified within 10% uncertainty and is consistent with other semi-automated methods [10], the BGS method identified 60% more individual pockmarks than GA1. This suggests that there are likely considerably more individual pockmarks in the overall study area, than the 220,508 pockmarks identified. Uncertainties associated with pockmark scour direction, although not quantified here, were considered low due to the large sample number (~40,000; supplementary material). However, for the 'slice' of directions associated with bearings having similar bearings to the survey lines (i.e., between 135 • and 150 • , and its opposite 315 • to 330 • ), this 'slice' may be misrepresented because step 4 of GA1 method eliminated these polygons ( Figure 5).
Considering that all three methods were developed separately and with different objectives in mind, the results are promising in terms of semi-automatically identifying pockmarks from a range of pockmark fields of varying density (Figure 7). The analysis suggests that while the objective of the study may drive the choice of a preferential method, e.g., scoured pockmarks may be better mapped using a method based on BPI algorithm (GA1), using a combined approach may be ideal to fully characterise any type of pockmarks. Moreover, the comparative analysis provides confidence in using the results of GA1 method over the whole study area to address hypotheses about pockmark formative processes.

Pockmark Formation and Ecosystem Contribution
Examination of shallow sub-bottom profiles across the survey areas suggests an absence of deeply-sourced gas with no evidence of gas migration from the sub-surface ( Figure 2). However, in the sediment samples examined, the concentration of CO 2 in sediment porewater (as TCO 2 ) is particularly high, approximately twice as high in these sediments compared to elsewhere on the Australian continental shelf (Figure 11). Without evidence for gas sourced from depth, this suggests that the CO 2 is derived from recent (i.e., shallow) sediments, and in turn must imply that organic matter (OM) trapped in the shallow sediment is breaking down rapidly. Our hypothesis for the origin of high pockmark densities in the Oceanic Shoals involves the concept of OM priming (or cometabolism), which pertains to the enhanced remineralisation of refractory sediment-bound OM through the breakdown of more labile OM (see conceptual model in Figure 15) [50]. Carbonate banks in the study area rise tens of meters above soft sediment plain where seabed sediments are often in the euphotic zone (~<60 m) and conditions thus are suitable for the growth of benthic algae as the local source of OM [37]. The shedding of labile organic matter from banks is a contemporary mechanism and likely occurs due to the relatively high current velocities in the region. This is supported by the hydrodynamic current model results (Figures 13 and 14) and the observed scouring of the pockmarks (Figure 3). The lower TOC concentrations per unit surface area of sediment in the high-density pockmark fields may provide evidence for OM priming, on the assumption that some particle bound OM had been mineralised into dissolved CO 2 (Figure 10b). In our conceptual model (Figure 15), the high porewater CO 2 concentrations observed reflect a combination of factors including elevated sedimentary TOC concentrations, the break-down of labile OM shed from the banks, the cometabolism of particle-bound OM and restricted permeability imposed by fine sediments (including the marl-like CaCO 3 content in the mud fraction). Pockmarks could arise if this combination of factors creates saturation with respect to CO 2 followed by gas production.
The data from the Oceanic Shoals region are similar to data from the Petrel Sub-basin ( Figure 1) where high TCO 2 in shallow sediments and pockmark formation was attributed to the breakdown of OM and CaCO 3 in the shallow subsurface [29,33]. There, pockmarks are abundant and their distribution correlates strongly with the known distribution of multiple shallowly buried paleochannels, similar to Area 2 here (Figure 2). In the Petrel Sub-basin study, mangrove-derived OM was abundant. Whilst in this study the OM is understood to be sourced from marine biota that inhabit the banks, the co-location and similarity in shape of the high-density field with the cross-shelf valley suggest that OM may also be sourced from the buried paleochannel. Although the OM source generally appears different for the study area (mainly banks), there is a broad similarity in pockmark formation in these two separate seabed regions of the Bonaparte Basin. This hypothesis, in turn, suggests that it likely can be applied elsewhere in this immediate region, including the pockmark fields identified on the Darwin Shelf [24]. It has been suggested that pockmark formation in this region may be due to the presence of banks focusing fluid flow from the shallow sub-surface [29]. The evidence from this study neither supports nor negates the latter hypothesis.
The biogeochemical and physical environment of the pockmark fields appear to drive some of the polychaetes biodiversity in this area. Infauna analysis showed distinct polychaete assemblages and higher richness associated with low density pockmark fields (<25/km). These fields occupy 15% of the study area, are located randomly ( Figure 6) and are associated with coarser grain-size (large proportion of shelly sand and gravel texture, as per high % CaCO 3 ) where TCO 2 concentrations are at their lowest ( Figure 15). In some environments, sediment grain size can regulate infaunal communities, with sand content increasing oxygen permeability and associated habitat availability beneath the seafloor surface [51]. However, sediment grain-size is unlikely to be the primary determinant of infaunal species distributions, with other factors such as hydrodynamics, nutrient availability, and biological interactions being key drivers [52]. This is supported by previous analyses of polychaete assemblages in northern Australia in which sediment grain-size was shown to be weakly related to polychaete assemblages and richness [53,54]. Alternatively, as per the intermediate disturbance hypothesis [15], polychaetes may be responding to seabed disturbances associated with active pockmark fields, with successional communities being more prevalent in pockmarks as disturbances such as gas expulsion occur ( Figure 15). For example, two groups of pockmarks in the Bay of Fundy were each associated with discrete megafaunal communities, one category supporting pre-equilibrium communities and the other supporting equilibrium communities [53]. Further work is needed on this relationship, including targeted sampling within and outside individual pockmarks. Figure 15. Conceptual model summarising the hypothesis about the formative processes of the pockmark fields in the study area. (a) Overview of the bank-plain system; (b) Low pockmark density scenario where more elevated sand content results in lower porosity compared to finer grain-sizes (shown here as clay and silt particles) and presumably higher permeability (interconnected pore space). Therefore, when OM is mineralised into DIC it does not build up to saturation levels for CO2. Polychaete biodiversity is relatively high because of the increase gas exchange (O2) and more stable conditions. (c) Medium pockmark density scenario where porosity is higher than (a). Permeability may become occluded within sediment micro-niches due to the occurrence of fine sediments (silt and clays) and the swelling of smectite. Limited pockmark formation occurs due to gas build up. Polychaete biodiversity is likely low due to the low O2 concentrations and increasing CO2 concentrations caused by restriction in permeability and the physical disruption to the sediments due to pockmark formation. (d) High pockmark density scenario in which sediments contain more silts and clays creating more micro-niches for gas formation. The gas is expulsed and the seafloor becomes "rugged". Clay particles have less sediment-bound OM due to OM priming. Sediments are mobilised from pockmarks due to strong near-seabed currents. Low polychaetes biodiversity is presumably due to physical disturbance.
The data from the Oceanic Shoals region are similar to data from the Petrel Sub-basin ( Figure 1) where high TCO2 in shallow sediments and pockmark formation was attributed to the breakdown of OM and CaCO3 in the shallow subsurface [29,33]. There, pockmarks are abundant and their distribution correlates strongly with the known distribution of multiple shallowly buried paleochannels, similar to Area 2 here (Figure 2). In the Petrel Sub-basin study, mangrove-derived OM was abundant. Whilst in this study the OM is understood to be sourced from marine biota that inhabit the banks, the co-location and similarity in shape of the high-density field with the cross-shelf valley suggest that OM may also be sourced from the buried paleochannel. Although the OM source generally appears different for the study area (mainly banks), there is a broad similarity in pockmark formation in these two separate seabed regions of the Bonaparte Basin. This hypothesis, in turn, suggests that it likely can be applied elsewhere in this immediate region, including the pockmark fields identified on the Darwin Shelf [24]. It has been suggested that pockmark formation in this region may be due to the presence of banks focusing fluid flow from the shallow sub-surface [29]. The evidence from this study neither supports nor negates the latter hypothesis. density scenario where more elevated sand content results in lower porosity compared to finer grain-sizes (shown here as clay and silt particles) and presumably higher permeability (interconnected pore space). Therefore, when OM is mineralised into DIC it does not build up to saturation levels for CO 2 . Polychaete biodiversity is relatively high because of the increase gas exchange (O 2 ) and more stable conditions. (c) Medium pockmark density scenario where porosity is higher than (a). Permeability may become occluded within sediment micro-niches due to the occurrence of fine sediments (silt and clays) and the swelling of smectite. Limited pockmark formation occurs due to gas build up. Polychaete biodiversity is likely low due to the low O 2 concentrations and increasing CO 2 concentrations caused by restriction in permeability and the physical disruption to the sediments due to pockmark formation. (d) High pockmark density scenario in which sediments contain more silts and clays creating more micro-niches for gas formation. The gas is expulsed and the seafloor becomes "rugged". Clay particles have less sediment-bound OM due to OM priming. Sediments are mobilised from pockmarks due to strong near-seabed currents. Low polychaetes biodiversity is presumably due to physical disturbance.

Pockmark Maintenance and Contribution to Hydrodynamic Knowledge
Correlation of the observed pockmark scours ( Figure 5) with the hydrodynamic model results ( Figure 14) suggest that at a regional scale, seabed currents are dominated by the SE-NW tidal flow. Results also indicate that small differences in the predominant modelled and observed directions between each mapped area exist (Figures 5 and 14). Given this strong regional agreement, the fine-scale variations observed within an area can give valuable insight into seabed and oceanographic processes at a scale finer than model resolution (Figures 3 and 7). For example, the scattering of the scour directions observed in Area 1 compared to the more dominant bi-directionality observed in Areas 2 and 3 ( Figure 5) suggest that the presence of larger and shallower banks generates more turbulence.
Though the modern seabed within the Bonaparte Basin is in places sediment-starved [21,28], turbulence within pockmarks enhances sediment suspension [11,12]. Given the strong tides here, the presence of open pockmarks within the organic-rich seabed sediments of the Bonaparte shelf suggests that they are active and that there is insufficient sedimentation to bury them faster than they form. Cores collected within the central Bonaparte Basin to the southeast of this study area suggest that post-LGM sediment thickness in the center of the basin is less than few meters [26]. This may explain the small size and depth of the pockmarks with smaller pockmarks typically forming in thinner sedimentary strata [55]. Furthermore, even though high turbidity in the water column was observed [20], the overall bed shear stress is enough to support erosion and mobilisation, while preventing settling of the sediments (Figure 14).
The wide spectrum of shapes exhibited by the pockmarks indicates that they provide an integrated record of the various stages of their formation (Figure 7). It is however, difficult to attribute whether scouring is associated with the formation of the pockmark at a specific time within a tidal cycle or that it reflects a long-term scouring process attributed to the integrated tidal regime. Some studies show that eccentricity of pockmarks is made at the time of the pockmark formation [18] and that it is unlikely that modern currents are reshaping inactive pockmarks, while others suggest that pockmark-scale turbulence plays a role in the excavation and scouring of pockmarks [11,12]. The intermixed variability in the shapes and scour directions observed over a small area here indicates that currents at the time of formation drive the eccentricity of the pockmark, but that scouring is ongoing thereafter.

Conclusions
On the northern Australian continental shelf, pockmark formation appears related to the decomposition of organic material which is richly preserved within the shallow sub-surface. Where carbonate banks are the primary point source for organic detritus, pockmarks become concentrated around the banks and are modified by tidal currents that are steered by the banks. As such, the morphometric data extracted here for large numbers of pockmarks provides a robust proxy to support modelling of local to regional hydrodynamics. Local variability in sediment texture and geochemistry around carbonate banks also reflects local organic input and appears to influence polychaetes biodiversity within sediments; but it also relates to seabed disturbance as driven by pockmark formation. In summary, there are multiple interactions that influence pockmark form and distribution. As more pockmarks fields are unveiled by the increased availability of high-resolution seabed data, semi-automated characterisation of pockmarks will become essential and better understanding of the seabed processes will be achieved.