Discovery of numerous pingos and comet-shaped depressions offshore southwestern Taiwan

High-resolution bathymetry collected with an autonomous underwater vehicle (AUV) along the flanks of three ridges of the accretionary prism offshore southwestern (SW) Taiwan revealed more than 650 elongated depressions in water depths ranging from 1155 to 1420 m. The depressions are between 12 and 129 m long, 5 to 70 m wide, and up 9 m deep at their center and shallowing downslope to about 1-m depth. Due to their shape in downslope cross section, they are termed comet-shaped depressions (CSD). The CSD occur in patches of more than 100 with densities of 53 to 98 CSD/km2. In addition, seven topographic mounds were mapped and interpreted as pingos, which remotely operate vehicle (ROV) observations and sampling show to be covered with authigenic carbonate. These features overlie areas where multichannel seismic reflection (MCS) profiles show bottom simulating reflectors (BSR) and dipping strata extending from below the BSR to near the seafloor. We consider comet-shaped depression, a new type of pockmark, forms on a sloping seafloor where fluids expulsion occurred. We also suggest that the two types of distinctive geomorphic features are attributed to fluid venting which occurs at different rates, with the mounds developing slowly over time, but the CSD forming in discrete events perhaps associated with large earthquakes.


Introduction
Fluid expulsion at the seafloor could release large volumes of methane, influence submarine slope stability, and even accelerate global warming (e.g., Berndt et al. 2014;Evans et al. 1996;Cochonat et al. 2002). Several distinctive types of seafloor morphological features are commonly associated with fluid escape and are widely distributed along continental margins, such as pockmarks (Hovland 2002), mud volcanoes (Dimitrov 2002), and submarine pingos ). The role fluid migration plays in generating these morphological features is increasingly gaining attention, in part because such morphologies are helpful for identifying and understanding potential geohazard issues. The advent of near seafloor survey tools, such as deep-towed cameras, autonomous underwater vehicles (AUV), and remotely operate vehicles (ROV), enable the visualization of deep-sea micromorphology, provide new evidence of relatively small-scale fluid seepage structures in modern marine sediments, and improve the understanding of the underlying processes (e.g., Caress et al. 2008;Paull et al. 2008, Paull et al. 2015Wynn et al. 2014).
Evidence for fluid seepage exists offshore southwestern (SW) Taiwan in the form of mud volcanoes, pockmarks, authigenic carbonate build-ups, and chemosynthetic biological communities containing bacteria mats and vesicomyid clams Huang et al. 2011). In addition, multichannel seismic reflection (MCS) profiles show strong and widely distributed bottom simulating reflectors (BSR) (Chi et al. 1998;Liu et al. 2006). The presence of BSR suggests that free gas and gas hydrates exist in the subsurface. The occurrence of elevated methane concentrations measured within near seafloor sediments and in bottom waters is commonly reported offshore SW Taiwan Yang et al. 2006), suggesting that methane is seeping out of underlying methane-rich zones. However, some methane anomalies in sediment and water column samples occur in locations where no distinctive seafloor features are observed by ship-borne geophysical surveys Chuang et al. 2006Chuang et al. , 2010. In addition, high-resolution surveys suitable for identifying small morphological features associated with seepage at accretionary wedge are rare. Thus, the distribution of seafloor seepage and their relation to fluid flow pathways remain unclear along accretionary wedges offshore of SW Taiwan.
In 2013 and 2017, AUV surveys (~1-m lateral resolution seafloor data) were conducted over three ridges: Good Weather Ridge (GWR), South Yung-An East Ridge (SYER), and Four-Way Closure Ridge (FWCR) in the accretionary prism offshore SW Taiwan (Figs. 1 and 2). Mini-ROV dives were also conducted to provide ground truths visual images of features identified in AUV surveys. MCS profile data are used to outline the underlying structure in the study area. We aim to (1) identify and delineate seepage-related seafloor features and (2) discuss the formation of these seafloor seepage features. A conceptual model is presented that attribute the comet-shaped morphology to seafloor seepage which occurs on a sloping seafloor.

Geological background
Taiwan is located at a convergent tectonic boundary between the Eurasian Plate and the Philippine Sea Plate (Ho 1986;Huang et al. 1997;Teng 1990). The area offshore SW Taiwan is an incipient collision zone where the accretionary wedge is encroaching on the passive China continental margin ( Fig. 1) (Liu et al. 1997;Liu et al. 2004). The accretionary wedge is characterized by a series of ridges formed by folds and thrusts (Liu et al. 1997;Liu et al. 2004;Lin et al. 2008). The accretionary wedge is considered to be a fluid expulsion prone region and is known to host pockmarks, mud volcanoes, fluid seeps, and authigenic carbonates ( Fig. 1) (Klaucke et al. 2016;Lin et al. 2006;Lin et al. 2009;Schnürle et al. 2011). Evans and Fischer (2012) suggest that fluid may migrate along dipping permeable strata within the fold-and-thrust structures. In addition, the steepness of the slopes on the flanks of these ridges (3°to 15°; Lin et al. 2008) and frequent powerful earthquakes in the arc-continent collision belt (e.g., Lin et al. 2015) make the flanks of these ridges susceptible to slope failure.
The uplift associated with the actively forming fold and thrust topography and erosion has enhanced the morphologic relief in the region, between Penghu and Kaoping (aka Gaoping) submarine canyons (Fig. 1). According to Yu and Hong (2006), submarine canyons in this study area started forming during the Pliocene. In addition, Penghu and Kaoping Canyons are considered to be tectonically controlled (Liu et al. 1993;Yu and Hong 2006;Chiang and Yu 2006). While the canyons are major pathways for sediment transport from the upper slope to the basin leaving sediment (Lin et al. 2008, Lin et al. 2014, the troughs between the accretionary ridges ( Fig. 1) are slowly accumulating sediment. , a chirp side-scan sonar with 110 kHz center frequency, and a 1-6-kHz chirp subbottom profiler. The AUV was navigated using a Doppler Velocity Logger and an inertial navigation system . The AUV was programmed to travel at 3 knots while maintaining an altitude of 50 m above the seafloor and with~150 m line spacing to provide overlapping multibeam swaths. Multibeam data and side-scan sonar images were processed to 1-m grid resolution bathymetry using MB-System (Caress and Chayes 1996;Caress et al. 2008). Chirp data with 10-cm vertical resolution were interpreted using the IHS Kingdom Software. Slope and aspect maps have been generated to reveal the small-scale morphologic features and to help quantify their shapes ( Fig. 3a; Table 1). Two survey lines in the GWR area, initially conducted in 2013, were repeated in 2017.

ROV observations and sampling
Several dives of MBARI's mini-ROV were conducted at seafloor sites in May 2017 to provide visual ground truth for the AUV surveys and to collect samples using a manipulator arm and push cores. In this study, we present results from mini-ROV dive 97 at the SYER site. In total, 3-h of video observation of the seafloor, two partly lithified rock or crust samples, one biological sample, and four push cores were collected on a mound in 1240-m water depth. Subsamples of the lithified rock or crust samples were analyzed by Accelerator Mass Spectrometry for 14 C content by Beta Analytic Radiocarbon Dating Laboratory, Miami, Florida.

Multichannel seismic profiles
MCS profiles were collected by the R/V Ocean Researcher I in 2017. All MCS data were recorded by a 108-channel, 1350-m long streamer. The seismic source used is a threeairgun array with a total volume of 505 cubic inches, firing at 25 m shot spacing. The seismic data processing flows include SEG-D data input, trace editing, defining geometry, bandpass filtering, amplitude compensation, normal moveout correction, CMP stacking, and time migration.
All the data were processed using the ProMAX seismic data processing software, and then input into the Kingdom Software for interpretation.

Regional bathymetry
Regional bathymetry (100-m grid modified from Liu et al. 2004) in this study area shows several distinctive elongate ridges in 900 and 2000-m water depths trending NE-SW to N-S between the Penghu and Kaoping Canyons (Fig. 1). Between these ridges are five flat slope basins that deepen to the SSE (Fig. 1). The flanks of these ridges are associated with scarps that slope at up to 15°and frequently form arcuate embayments that range from~100 m to more than a kilometer across (Fig. 2).

Comet-shaped depressions
The high-resolution bathymetry (1-m grid) covering portions of the GWR, SYER, and FWCR ( Fig. 1) shows numerous depressions on their lower flanks in water depths from 1155 to 1420 m (Fig. 2). These depressions range in length from 30 to 70 m, are up to 9 m deep at their center, and become progressively shallower downslope. In cross section, these depressions have an asymmetric shape with a long tail trending downslope, producing a morphology that resembles a comet. We will use the descriptive name of comet-shaped depressions (CSD) for the remainder of this text (Fig. 3a).
On the west flank of GWR ( Fig. 2a), more than 150 CSD occur on the otherwise smoothly sloping lower flanks in 1155to 1240-m water depths between the steep upper scarp and the basin floor to the west. One patch of rounded depressions without a long tail is observed on the comparatively flat seafloor in the northeastern side of this AUV map (Fig. 2a). On both sides of SYER and FWCR (Fig. 2b, c), there are widely distributed CSD in water depths of 1250 to 1320 m and 1350 to 1420 m, respectively.
The morphologies of the CSD are similar within each patch and between the three ridges ( Table 1). The CSD occur where the surrounding seafloor slope is between 7°to 10° (Fig. 2). CSD are not observed on the more-gentle slope on the flat top ridges or within the sediment filled basins. The CSD are also restricted to areas where there is some sediment drape onlapping onto the flanks of the adjacent scarps.
The CSD occur as single entities, characteristically spaced 25 to 160 m apart from each other. Within the patches, their frequency of occurrence ranges from 58/km 2 in FWCR to 98/ km 2 in SYER (Table 1). Some small groups of CSD that seem to be aligned in rows occur where the truncated edges of dipping bedding planes occur in the near subsurface (Fig.  3b, c).
The morphologic expression of all the CSD is consistent within and between the three areas. No aureoles or rims are observed that might suggest deposits of excavated material on their flanks, but subtle tails extend directly down slope below some CSD (Fig. 3b, c). In particular, side-scan sonar image reveals similar backscatter intensity within the CSD patches (Fig. 4b), suggesting that they are all at similar stages of morphological development.

Circular topographic mounds
The high-resolution bathymetry covering portions of the GWR, SYER, and FWCR ridges (Fig. 1) also show seven circular topographic mounds with distinctive rough surface textures (Fig. 4). On the west side of the GWR, a mound of 25 m in diameter and 14-m relief height is present in 1212-m water depth near the CSD zone (Fig. 4a). Two mounds about 100 m in diameter in the SYER are not only revealed by highresolution seafloor morphology but also show high backscatter intensity on side-scan sonar images (Fig. 4b-d). A large mound which is 140 m in width, 350 m in length, and 10 m in height in 1340-m water depth is observed on top of the FWCR (Fig. 4e). Side-scan sonar images show water column reflections above the rim of the large mound, indicative of rising gas bubbles. Three similar mounds are identified on the northwestern side of the FWCR near where CSD are also observed (Fig. 4f). Two of these mounds occur within or on the flanks of CSD depressions (Fig. 4f).
ROV observations show that the mound in the SYER is largely covered with authigenic carbonate which occurred as large blocks (Fig. 5a, b). The carbonates appear to have originally formed as a pavement within the typically seafloor sediment, which has now broken into the blocks. The carbonates contained cemented clam shells which are inferred to be vesicomyids (Fig. 5c). Pore water are extracted from push cores contained elevated sulfide concentrations (S. Lin, personal communication). Some patches of tubeworms (inferred to be paraescarpia) are also seen and sampled (Fig. 5d) (S. C. Chen, personal communication). Both the clams and tubeworms are characteristic for methane seeps (Kiel 2010).
Measurements of the 14 C content in two samples of the carbonate collected from the mound yielded concentrations of 1.45 ± 0.04% modern carbon (pMC) and less than 0.44 pMC, which are close to or beyond the detection limit, respectively. The δ 13 C values of the same samples are − 39.9 and − 44.3 per mil VPDB, respectively, showing substantial proportions of methane-derived carbon.

Sub-bottom characteristics
Both the chirp profiles and MCS profiles show that the CSD occur where the drape of sediments from the trough thin and lap onto the flanks of the ridges (Fig. 6). The chirp profiles show little or no sediment drape (< 30 cm) has accumulated since the CSD formed (Fig. 6a). Moreover, the chirp profiles do not reveal any older, now buried CSD in the shallow strata.
Dipping parallel reflections are seen in MCS profiles that project to near the seafloor where CSD occur (Fig. 6). The reflections that extend across the BSR pass through zones with high amplitude reflector packages (HARPs) and terminate near the seafloor where these CSD occur (Figs. 6b). These dipping reflections are truncated on both flanks of the anticline. Some small subsurface fractures occur beneath the CSD on the eastern flank of the SYER (Fig. 6b). In addition, the MCS profiles show that a chimney structure is present beneath the mound and extends upward from below the BSR (Fig. 6b).
Subtle differences in two generally similar MCS profiles crossing the flank of the GWR (Fig. 7) help explain the distribution of CSD. Slightly dipping HARPs with reversed polarity extend to the seafloor where CSD occur (Fig. 7a). However, the MCS reflection profile (Fig. 7b) shows that dipping HARPs with reversed polarity do not extend to the seafloor along the flank of the ridge, but instead to where circular depressions were mapped on gently sloping seafloor on the east side of the prominent ridge (Fig. 2a).

Discussion
The formation of seafloor seepage features

Circular topographic mounds
The seven mounds identified in the bathymetry are inferred to be related to seepage. The large chimney structure (Fig. 6b) suggests a sub-surface pathway for fluid movement. The existence of chemosynthetic megafauna from the one sampled mound shows that seepage is on-going (Fig. 5). Flares have been reported over other mounds which indicates active seepage (Klaucke et al. 2016). The low δ 13 C values and largely depleted 14 C content of the authigenic carbonate confirm the presence of methane in the seeping fluid and suggest that it has been occurring for some time.
The occurrence of the distinctly elevated topography associated with these mounds needs to be explained. These features do not show either the flow textures or presence of extruded material characteristic of mud volcanoes visited on other mini-ROV dives off Taiwan (Hsu et al. 2017). Moreover, the pavements of authigenic carbonate indicate a lateral continuity of the seafloor, which is inconsistent with these features being erosional remnants or accreted chemoherms (e.g., Teichert et al. 2005). However, these mounds are geomorphically similar with mounds from methane seepage site elsewhere which have been interpreted to be gas hydrate pingos (Hovland and Svensen 2006;Paull et al. 2008Paull et al. , 2015. Gas hydrate pingos are mostly reaching a height from 10 to 40 m above the seafloor and 100 m or more in width (e.g., Paull et al. 2007;Waage et al. 2019). The relief on the mounds we observed are in water depths of 1212-1340 m reaching a height of up to 10 m above the seafloor and are up to 350 m wide. Thus, these mounds are interpreted to be associated with gas hydrate growth in the near subsurface which has blistered the seafloor.

Comet-shaped depressions
Widely distributed CSD are recognized on the flanks of these ridges from ultra-high-resolution bathymetry (Figs. 2 and 3). To our knowledge, such features have not been reported previously in gas hydrate system, perhaps because they are too small to be detected by conventional ship mounted geophysical systems on such high slopes. However, several mechanisms have been proposed for the formation of seafloor   depressions including (1) erosion by bottom currents (e.g., Hillman et al. 2018;Klaucke et al. 2018;Loncke et al. 2016), (2) slope failures which leave arcuate landslide scars (e.g., Baeten et al. 2013;Hampton et al. 1996), (3) spring sapping which recharges the platform fluids by entraining seawater (e.g., Paull et al. 1990, Paull et al. 1991, and (4) fluid venting which typically leaves circular depressions commonly called pockmarks (e.g., Harrington 1985;Hovland 2002;Judd and Hovland 2007). No evidence exists for the existence of strong bottom currents in the areas directly overlying the CSD. Moreover, bottom current at these water depths would flow along the contours. However, the elongation of the CSD consistently trends downslope (Fig. 2b, c), which is perpendicular to any inferred contour currents flows. Thus, bottom current is unlikely to be a major factor in the formation of the CSD.
Another possible mechanism to cause similar seafloor depressions is landslide (Hampton et al. 1996). Several large scarps occur on the flanks of three ridges (Fig. 2) which have been interpreted by Klaucke et al. (2016) to be landslide scars, but the size of the CSD is too small to be identified in his data. The CSD could be smaller-scale landslides (Baeten et al. 2013). However, why the CDS occur in fields, why they occur as single entities which are characteristically spaced 25 to 160 m within the fields, why they are so similar to each other, and where the material went are unclear and not explained by there being simply small landslides. Spring sapping also could provide comet-like features. However, the Florida Escarpment (e.g., Paull et al. 1991, Paull et al. 2002 which is erosion by spring sapping is limited as the primary mechanisms by chemical dissolution. Besides, spring sapping is effective on limestone; it would be ineffective within the hemipelagic and clastic sediments in the study area. The potential role of fluid venting on the formation of CSD should be considered. Fields of pockmarks (e.g., Harrington 1985;Hovland 2002;Paull et al. 2002;Judd and Hovland 2007) characteristically occur in areas where the seafloor slope is relatively low (e.g., < 1°), and there is a thick Holocene sediment drape. The main morphological difference between pockmarks and CSD is the long tail of CSD which are oriented directly downslope. If material was excavated to form these depressions by fluid venting, it presumably drained down slope as a suspended sediment cloud, leaving no significant trail of intact debris. The areas where fields of CSD occur are coincident with the occurrence of dipping strata form fluid flow conduits connecting into the subsurface and relatively steep slopes (7°to 10°; Fig. 2). We infer that the CSD is a new type of pockmark that forms on a sloping seafloor (Fig. 8) instead of the more common rounded shape pockmarks.  Fig. 7b is indicated by SSS label. Gray lines indicate some fractures on the profile Fig. 8 A conceptual sketch of pockmarks and comet-shaped depressions (CSD), outlining the suspected fluid pathways. Fluid accumulation beneath the seafloor and blowout vertically from zones of weakness to form the pockmark in the flat region. Instead, rising fluids along higher permeability layers which are dipping strata on the flanks of the ridge and vent on a slope to form the CSD Fig. 9 A conceptual sketch of small-scale seep evolution in anticline and monocline. a Initial, structural uplift, and landslides occur. b Before CSD formed, some dipping high-permeability layers exposed near a slope seafloor and thin sediment draped. c CSD formed, rising fluids vent on a slope seafloor to form the CSD during the variation of hydrostatic pressure event. Dark gray layers indicate the low-permeability layer, and light gray layers indicate the high-permeability layer. Cyan arrows show the pathway of fluid along a dipping high-permeability layer. Cyan area on the top of the ridge indicates the fluid accumulation What are the main factors to form CSD in accretionary ridges?
The fluids seeping onto the seafloor in accretionary wedges commonly migrate upwards along permeable conduits like thrust faults (e.g., Evans and Fischer 2012;Judd and Hovland 2007;Barnes et al. 2010;Klaucke et al. 2016). The dipping HARPs extend to near the seafloor where CSD field occur (Figs. 6b and 7a) and are likely associated with high permeability strata, containing gaseous methane, and serve as conduits (e.g., Benjamin and Huuse 2017;Wenau and Spiess 2018) which would direct fluids to the areas where CSD occur.
We infer that the existence of the underlying dipping and high permeability layers which provide a fluid pathway is an important factor in the formation of the CSD (e.g., León et al. 2014). This interpretation explains why CSD are not observed on the ridge tops and one of the northwest-dipping smooth slopes in GWR, as high permeability layers do not extend close to the seafloor there ( Figs. 2 and 7). That also explains why some CSD appear to be aligned along the truncated edge of shallowly buried bedding planes (Fig. 3b, c). However, the chirp profiles within the CSD fields do not show acoustic blanking, a common feature of active seepage zones (e.g., Papatheodorou et al. 1993;Zitter et al. 2008). Thus, there is no evidence that gaseous methane is present within the sediment drape at this time.
When did the gas hydrate pingos and CSD form?
While the available data do not define the formation time accurately, we suggest that the inferred gas hydrate cored pingos are still actively growing. The observation of active chemosynthetic communities and rising gas bubbles in the large mound of the FWCR indicates that methane-bearing fluids are still seeping to the seafloor. The growth of up to 14-m high pingos mounds and the occurrence of large blocks of slowly forming 14 C-depleted authigenic carbonate containing fossilized clams suggest that sustained venting has occurred at these sites (Aloisi et al. 2004;Luff et al. 2005) (Figs. 4 and 5). Therefore, the pingo mounds may require thousands if not tens of thousands of years for their formation.
In contrast, we infer that these CSD could be intermittently active, which appear to have formed sometime in the Holocene. Firstly, these CSD do not appear to be active because no flare signals have been detected above them and there is no acoustic blanking within the adjacent sediment drapes. In addition, even though there are twelve earthquakes between magnitude 2.9 and 4.5 occurred between 2013 and 2017, no differences were noted in the bathymetry in the small area of overlap during the 5-year interval. Secondly, the rather consistent 4-m depth of the CSD (Table 1) and absence of a significant sediment drape cover provides some constraint on their age. Using the lowest sedimentation rate estimate for this area (0.7 m/kyr; Lin et al. 2014), the maximum age of these sediment into which the CSD cut is~5.7 kyr. However, the lack of significant sediment fill suggests a considerably younger age.

Evolution of CSD
Based on all the observations and data analyses, we propose a simplified conceptual model to explain the CSD formation on the sloping flanks of accretionary ridges. Structural deformation and uplift within the accretionary wedge produce anticlinal and monoclinal ridges (Figs. 6b and 7a). The steepening flanks of these ridges become susceptible for failure. Periodic large landslides formed the arcuate scarps on the steep sides of the developing anticlines and monoclines and exposed the underlying dipping and high permeability layers (Fig. 9a). With time, hemipelagic sediment gradually accumulates on the slide scars covering the exposed high permeability conduits with a drape lapping onto the lower flank of these ridges and progressively thin upslope (Fig. 9b). The CSD are simply pockmarks which developed within the sediment drape on these steep slopes when fluids vented as a suspended sediment cloud through the permeable conduits.
The similarity of morphology between CSD within and between the CSD fields on the three ridges, the limited sediment fill, and absence of older or buried CSD suggest that they were all formed within a narrow range in time, indistinguishable from being associated with one discrete region-wide event. However, we cannot eliminate the CSD being a composite of processes which occurred over a century or so. We cannot rule out that the CSD can be recurrent or are still developing. In other words, no evidence of earlier cycles suggests an explanation of what those cycles would be. Therefore, given the limited sediment cover, this would suggest that the latest event occurred sometime in the late Holocene.
The event presumably involved enhanced flow, perhaps even vigorous flow, through the underlying permeable conduits to excavate the depressions. Such an event would be driven by region-wide pore pressure variations as it influenced all three ridges in a similar fashion. Strong compressive stress, which occurs during large earthquakes (e.g., Lin et al. 2015), is a mechanism that might produce significant pressure transients which could cause significant flow through the dipping horizons (Fig. 9c). Such a flow is inferred to disturb the poorly consolidated near seafloor sediments which draped over the truncated edges of these dipping permeable horizons.
Despite Taiwanese historical records that document earthquakes less than 200 years, if an extreme shock event is a trigger, some historical tsunamis induced by earthquakes could be considered, such as the disastrous tsunami that occurred in Kaohsiung and Pingtung in southern Taiwan at 1781 AD (Lau et al. 2010). Nevertheless, hitherto, no direct evidence of dating data collecting in the CSD. Therefore, we expect more observation and dating data for further analyses to understand what those cycles of intermittently active events would be.

Conclusions
High-resolution surveys reveal two distinctive types of geomorphic features on the flanks of three ridges in the accretionary wedge offshore SW Taiwan. Both are attributed to fluid venting. One type consists of seven authigenic carbonate covered mounds interpreted to be gas hydrate cored pingos, features observed previously in other methane venting environments. The other type is represented by more than 650 cometshaped depressions (CSD) that exist in water depths ranging from 1155 to 1420 m and with 50-m average lengths, 24-m average width, and 4-m average relief. They occur on steeply sloping ridge flanks with sediment drapes. We infer that CSD are a type of non-conventional pockmark which could be a widespread phenomenon in other accretionary ridges where abundant gases or gas hydrate presents with some high permeability layers and dipping fluid pathways exist beneath. Our data also imply that fields of CSD formed sometime in the Holocene, potentially in a single event. We suggest that the fields of CSD may be formed on a sloping seafloor by dewatering events which occur during powerful earthquakes.