Gas Emissions in a Transtensile Regime Along the Western Slope of the Mid-Okinawa Trough

Gas emissions from the seabed are favored by tectonically active settings and their distribution is often linked to the nearby faults. Here we use the multi-beam echo-sounder (MBES) and the multi-channel seismic (MCS) data and a sediment core to show multiple gas emissions near the fault complex out of the shelf of the Mid-Okinawa Trough. The features indicating the gas emissions include 1) a set of the conical positive reliefs at the seabed, 2) the bundle-shaped clusters of the high-backscattering intensities in the water column, and 3) the sub-circular medium-to high-backscattering patches at the level of the seabed. These features together show that the free gases can escape from the marine sediments then rise in the water column at present, while some other gases trapped in the sub-seafloor sediment might contribute to the precipitation of the authigenic carbonates in the past. The spatial relationship between the gas emissions and the faults suggests that the faulting driven by the back-arc extension should provide the permeable migration pathways for the gas emissions to operate, and thus determines where most of them could potentially occur. The area surrounding the restraining bend concentrates part of the gas emissions rather than along the fault lines, due to the lateral compression and the structural complexity. This is demonstrated by the results of the numerical model of finite element method (FEM), which shows two gas emissions are within the compressed zone of the modeled restraining step-over. This study provides new evidence of the role of the tectonic stresses in determining the sites of degassing of marine sediments.


INTRODUCTION
Gas emission in the marine environment has attracted significan research attention over the past few decades. The abundance of gas emissions could indicate the potential accumulation of hydrocarbons and, if ambient conditions were proper, the gas hydrates (Riedel et al., 2018). Gases leaked from the seabed constitute an important part of methane transported into the ocean (Judd et al., 2002;Feng et al., 2018). As a member of the potent greenhouse gases, methane can sometimes reach the shallow seawater by forming ascending streams of gas bubbles and possibly enter the atmosphere (McGinnis et al., 2006;Westbrook et al., 2009). This fate of methane has been questioned by the study in the Arctic Ocean, in which enhanced methane concentration was detected close to the gas emissions near the seabed but not in the seawater near the sea level (Myhre et al., 2016). Therefore, whether methane escaping from marine sediments can contribute to climatic warming is controversial. If dissolving into the seawater is the fate for most of this methane, the oceans would likely experience deoxygenation and acidification (Biastoch et al., 2011;Yamamoto et al., 2014).
Gas emissions emerge intermittently surrounding the continental shelves and slopes, and are sometimes coated by hydrates in the deep-water settings (Somoza et al., 2003;Sauter et al., 2006;Römer et al., 2012;Smith et al., 2014;Rümer et al., 2017). Free gases can bypass the less permeable sediments alone or be expulsed together with pore fluids and mud-rich sediments (Gennari et al., 2013;Cartwright and Santamarina, 2015). The geological records of such processes include a spectrum of surface expressions such as pockmark (Hustoft et al., 2009;Wang et al., 2018), mud cone (Roberts, 2001), mud volcano (Huguen et al., 2004) and carbonate mound (Schmidt et al., 2005). Gas emissions through these morphological features sometimes have the geophysical signature of the acoustic flares in the water column (i.e., gas plume, Heeschen et al., 2003;Dupre et al., 2015), which can be detected by multi-beam echo-sounder, and gas chimneys imaged by seismic data in the subsurface (Cartwright and Santamarina, 2015;Li et al., 2017).
Ongoing research probes into how geological processes control gas emissions and the one that has been studied most on global scale is the tectono-structural process (Ciotoli et al., 2020). Gas emissions have been detected in the extensional (Plaza-Faverola and Keiding, 2019), the compressive (Reed et al., 1990;Watson et al., 2019), the sheared (Huguen et al., 2004;Geli et al., 2008) and the mixed deformation regimes (Plaza-Faverola et al., 2014). The faults formed in these settings can indicate the stress regime surrounding gas emissions, and sometimes served as the escaping conduits for them (Ciotoli et al., 2020). The gas emissions in the Mid-Okinawa Trough have been inferred only by geochemical analysis (Xu et al., 2018) and what geological process controlled their distribution is not fully understood.
In this study, we use multi-beam echo-sounder (MBES) and multi-channel seismic (MCS) data to show multiple gas emissions out of the shelf of the Mid-Okinawa Trough and consider how faulting determines their potential locations. Here the nearby fault complex is formed in the transtensile regime, accompanied by local compressions surrounding the restraining bend. The finite element method (FEM) is used to simulate the stress regime surrounding two of these emissions in a restraining step-over. The consistency of the modeling result with the geophysical observation supports the role of local lateral compression in promoting gas emissions.

GEOLOGICAL SETTING
The Okinawa Trough is an incipient back-arc basin in the East China Sea and displays as an arcuate geometry in plan view ( Figure 1). It formed after the middle Miocene due to the crustal stretching (Lee et al., 1980), driven by the subduction of the Philippine Sea Plate under the Eurasian Plate (Kimura, 1985; FIGURE 1 | Bathymetric map showing the study area in the Okinawa Trough. The bathymetric data are downloaded from https://maps.ngdc.noaa.gov/viewers/ bathymetry (color bar-GEBCO2014 Bathymetry Color Scale). The methane seepages have been recorded at the sites of Cores A and C, C01, C10 and GT-D1 (Li et al., 2015;Sun et al., 2015;Xu et al., 2018).
Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 557634 Gungor et al., 2012). This tectonic activity accounts for the multiple periods of rifting and the latest one is still active at present (Kimura, 1985;Letouzey and Kimura, 1985;Yamaji, 2003). The en-echelon grabens bounded by the faults and the stepped, sometimes unclear rifting centers in the Mid-Okinawa Trough have been revealed before by some sparse seismic profiles (Kimura, 1985;Letouzey and Kimura, 1985;Sibuet et al., 1987). Such shallow structures and basaltic ridges have also been seen in the area of volcanic arc-rift migration phenomenon (VAMP, shown in Figure 1) (Sibuet et al., 1987). To the north, the faults formed by transtension have been mapped in detail by combining the multi-beam echo-sounder and the multi-channel seismic data in an area of ∼5300 km 2 (Figures 2A, 3). Their strikes change from ∼N60°E in the shelf-slope setting (marked by n1) to ∼ N70°E in the deep-water basin (marked by n2, Figure 2). The faults in the deep-water basin are mostly normal, while those in the slope often have scarps that are subparallel with each other (Figures 2,  3). There are some small-scale faults affected by strike-slip in the fault blocks, leaving minor positive relief at the seabed ( Figures  3F,G). The observations of the fault pattern are consistent with the previous interpretation of the basin-scale tectonic regime that is characterized by extension near the rifting center and strike slip in outer shelf-slope setting ( Figure 2B) Liu et al., 2016). The sedimentary study in the Mid-Okinawa Trough focuses on the sediment provenance and its variation during late Quaternary (Dou et al., 2010a, Dou et al., 2010bLi et al., 2019). The largest contribution to the sediment deposits there is from the Yellow River since the last deglaciation, while other limited sources include Yangtze River-derived and Taiwanderived sediments (Dou et al., 2010b;Li et al., 2019). In the Middle to Late Holocene, the sediments delivered into the basin was hindered by the barrier effect of the strengthened Kuroshio Current (Dou et al., 2010b;Li et al., 2019). The geophysical observation associated with sedimentary features in the Mid-Okinawa Trough is the shelf-margin delta clinoforms straddling around the shelf edge formed during sea-level lowstands in the shelf-slope setting (Berné et al., 2002;Li et al., 2014;Li et al., 2020). No considerable canyons and basin-floor fan have been detected yet in the area between 28°30′N and 29°10′N (Li et al., 2020). The near-seafloor sediments are mostly mud and the clay mineral is dominated by illite and smectite (core Oki02, Zheng et al., 2014;core M063-05, Li et al., 2019). An area of ∼1.2 × ∼0.3 km crust of the authigenic carbonate at the water depth of ∼600 m has been found near seabed methane seepage by the visual data recorded by a remotely operated vehicle (ROV) ∼200 km to the north of the study area. These outcropping carbonates consist of aragonite, pyrite and gypsum, and the isotopic measurements show that the methane feeding their growth is mostly biogenic (Cao et al., 2020).
The Okinawa Tough is known for the seafloor methane seepages and multiple ones along its western slope have been revealed by the record of the geochemical data from the pore water of the gravity cores (C01, C10, cores A and C, and GT-D1, marked in Figure 1) (Li et al., 2015;Sun et al., 2015;Xu et al., 2018). The relatively thin Sulfate-Methane Transition Zone at some sites (e.g., C25 in Figure 1) suggests rapid transport of dissolved methane toward the surface (Xu et al., 2018). Escaping of free methane gas has been inferred before to occur near the suspected mud volcano, while the geophysical evidence is recorded in the seismic cross section ( Figure 3). The features associated with free gases include the ∼2.2 km-wide acoustic turbidity out of the study area and some local enhanced reflections at different horizons within it ( Figure 3).

Acoustic Data and Sediment Core
We installed a Kongsberg EM122 multi-beam echo-sounder (MBES) system for the acquisition of the acoustic data of the seawater and the seabed. The system mounted by R/V Haidahao has a swath angle of 130°and can transmit and receive up to 288 simultaneous beams with a beamwidth of 1°. It has a sounding frequency of 12 kHz and up to 864 soundings per ping in the dualswath mode. The positioning and the monitoring systems are NAVCOMSF-3050 and Kongsberg SIS, respectively. The parameters of roll, pitch and heading have been corrected during a Sea Acceptance Test before data acquisition. The distances between the dip-oriented and the strike-oriented track lines are 4 and 8 km, respectively. The recorded data, vessel information, and sound velocity profile (SVP) have been integrated and processed by using the QPS Qimera software (including sound velocity correction, creating dynamic surface, etc.). The velocity of the seawater was measured by using the system of the Sea-Bird 911 plus CTD and the X-CTD. QPS Qimera and Fledermaus softwares were used for visualization of the gas plumes in fan view and space, respectively. The swath width increases from 0.94 km on the shelf (water depth of 200 m) to 5.9 km in the deep-water basin (water depth of 1090 m). The recorded backscatter values of the seafloor mosaic images were normalized in Fledermaus software. The cell size of the backscattering images is 11.83 m.
A SIG pulse L5 sparker (20-1000 Hz) was employed at the 1-2 m below the sea level for acquiring the two-dimensional (2-D) high-resolution seismic data. The ship towed the 48-channel streamer (294 m long). The data were sampled at an interval of 0.5 ms and then processed on board by using the ProMAX system (including pre-filtering, NMO, de-noising, etc.). The 2-D multichannel seismic (MCS) data can show the clear image of ∼550 ms in TWTT below the seabed with a vertical resolution of ∼2 m. Positive polarity is defined by a peak on the seismic trace and displayed as a black-red loop on the cross section.
The acquisition of the in-situ sediment cores at D5 was carried out by the seabed drilling rig Seabull. It was deployed from R/V Haidahao to recover sediments as deep as 55.1 m below the seafloor. The lithology of the sediment cores is not fully recorded.

Finite Element Method
The FEM has been widely used in the tectonic analysis by considering stress transferring across the interface of the neighboring elements and solving the matrix equation. In this study, we use the FEM to show the role of the shear stress in causing the local compression surrounding the restraining bend between two faults. The 2-D square plane is used to model the marine sediments in plan view, through which the  Figure 9. The area marked in grey color is set to eliminate the boundary effect. Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 557634 5 gas-rich focused fluid flow passed. The dimension of the model is 0.1% of the real one observed in the geophysical data (shown in Figure 4; Table 1) and the side length is 22.58 m. The peripheral 3.5 m-wide area is modeled to eliminate the effect of the boundary conditions (Figure 4). The sediments are assumed to be isotropic. The square plane was discretized into 18620 triangular elements (9508 nodes), most of which are equilateral ones (the discretization of the code distmesh2d can be downloaded from http://persson.berkeley.edu/distmesh/) and the side length is 0.25 m (Figure 4). The stiffness matrix of each triangular element [k], for unit thickness, is: where [B] is the strain matrix and [D] is the elastic matrix (Timoshenko and Goodier, 1982). They are expressed as: where E and υ are Young's modulus and Poisson's ratio, respectively.
A is the area of each triangular element and expressed as: Based on the principle of minimum potential energy, the global balance equation for the entire planar sediments is: where [K] is the global stiffness matrix and can be assembled from the stiffness matrix of each triangular element. u is the nodal displacement vector and in this study has two dimensions u 1 and u 2 (x-and y-direction, respectively). P is the nodal force vector. The displacements u 1 at left and right boundaries and u 2 at upper and lower ones are set to zero (Figure 4). This boundary condition makes the solutions of {u} are unique. Either of the fault zones is 8.88 m long and the horizontal separation is 1.58 m. The geometry of part of the fault is not input into the model so that the dense elements can reveal the features within the bend. Two groups of initial forces (F i ) are set to act on the nodes of two faulted zones and either group has opposite directions ( Figure 4). Equation 6 is solved using MATLAB's MLDIVIDE function. The nodal stresses are the weighted average of the neighboring triangular elements. The maximum and minimum principal stresses (σ 1 , σ 3 ) and the Von Mises stresses (σ VM ) are expressed as: where σ x , σ y and τ xy are the normal stresses along the x-and y-directions and the shear stresses, respectively. Similar to previous study (e.g., Liu et al., 2010), this model has been tested in a simple case of stepover (90°neutral). In this study, the modeling results are rotated counterclockwise by 39°for comparison with the geophysical observation. The modeling results of the absolute values of stress and displacements are uncertain but instead show a relative pattern. The physical properties of the rocks in the study area have not been tested before, and their values used here are from the previous study using FEM to model stresses in a general case (

Morphological Features
The seabed in the study area dips southeastwards and the dipping angle decreases seawards, ranging from ∼1.5°at the water depth of ∼700 m to ∼0.5°at >950 m. It is marked by a set of the aligned breaks in the relief along the fault scarps ( Figure 5). The linear-curvilinear scarps have a NE-SW trend and their separation distance in plan view is ∼2.5 km.
The vertical displacements of the scarps are 15-90 m and the offsets along each fault plane can be clearly observed on the seismic profiles (Figure 3). Some positive reliefs (named with Db) are located close to two of these faults (Figures 5, 6, 7).  Only one positive relief (Db29) is observed to the east of the nearby fault scarp (Figures 5, 6). The water depth, the length of the long axis and the elevation of the positive reliefs are shown in Table 2. Most of them have conical geometries ( Figures 5C,D). There is a ∼1 km long ridge-like positive relief elongated along the fault (Figures 5, 6). In the northeastern part of the study area (white box in Figure 5A), part of the NE-SW oriented fault scarp bends landwards. Db21 and Db29 are located at either side of the bending zone and have a similar distance away from it, forming a sub-symmetric feature in plan view (Figure 7).

Hydroacoustic Features
Five bundle-like clusters of the acoustic anomalies have been identified in the study area ( Figure 6). The gas plumes (named with Pb and numbers) have a higher backscattering intensity than the seawater of the background (Figure 6). The relatively highbackscatter intensities occur in the lower and the middle parts of each gas plume ( Figure 6). The gas plumes have a height of 163-355 m. Two of the gas plumes (Pb5 and Pb49) emanated from the surrounding areas of the positive reliefs, while the other three ones (Pb1-3) cluster around the arcuate part of one fault (Figures 5, 6).

Backscattering Features at the Seabed
There are some medium-to high-backscatter areas at the level of the seabed (shown as the bright black patches in Figures 5B,  6B). Their intensities are 5-10 dB higher than the mean of those of the background marine sediments. The bright patches are sub-circular, sub-elliptical or irregular in plan view ( Figure 5B). Most of them cover an area of <0.1 km 2 and the largest one occurs near the ridge-like positive relief (∼0.3 km 2 , Figure 5B). Part of the bright areas are located close to the fault scarps, while the rest occur at the places where there is no topographic anomaly ( Figures 5B, 6B). Not all of the positive reliefs are superimposed by the areas of the enhanced backscattering intensities.

Sediment Cores
D5 is located within a medium-to high-backscatter area of the seabed (Figure 8). The upper 6 m-long sediment cores at D5 mainly consist of greyish clay and a strong H 2 S odor was smelt during processing the samples prior to the on-board test. Angular shells occur at the depth of ∼0.45 m (Figure 8).
Neither gas hydrates nor soupy structures were found. Carbonates occur at the depth ranges of 1.33-1.85 m and 2.50-3.05 m deep and hollows and shells are common within them (Figure 8).

Numerical Model
The variations of the displacements along x-and y-directions (u 1 , u 2 ) display different symmetric patterns (Figure 9). The u 1 values are near zero in the middle zone, increasing and decreasing along the positive and the negative y-directions toward the upper and the lower boundaries. There are some anomalous values near the fault tips ( Figure 9). The variations of the u 2 values are symmetric about the faults and their linkage zone. Positive and negative values occur to the left and the right side of the symmetric line ( Figure 9). The normal stresses along x-and y-direction (σ x , σ y ) show the opposite variation pattern (Figure 9). The lobe-like zones of the positive and negative σ x values are located inside and outside the linkage zone, respectively. This pattern is converse for the σ y values ( Figure 9). The zones of high values of maximum principal stresses (high-σ 1 zone) and low values of minimum principal stresses (low-σ 3 zone) are both near the fault tips. For σ 1 and σ 3 , two zones extending from the fault tips coalesce within the linkage zone to form a sigmoidal geometry (Figure 9). The extending directions are slightly different between σ 1 and σ 3 . The σ 1 directions between the faults are almost perpendicular to that of the segment between them (red dashed line in Figure 9). The arrays sharing such direction together form a lozenge-like zone (blue dashed line in Figure 9). The zone of similar geometry also occurred in the modeling results of a transpression zone in the strike-slip setting by Nabavi et al., 2017. The directions of the maximal principal stresses within it will rotate more parallel with the trend of the fault as the underlap distance decreases (Nabavi et al., 2017). The lobes of enhanced Von Mises stresses (σ VM ) do not deflect near the fault tips and the coalesced zone shows a symmetric pattern about the segment linking the fault tips.

Acoustic Features Associated with Gas Emission
The gas plumes are compelling evidence of the present-day emission of the gas bubbles from the seabed. Their presence indicates that the free gas bubbles can rise through the water column by hundreds of meters before breaking up. It has to be stated that the frequency of 12 kHz used in this study for data acquisition is not high enough to reveal all of the emissions of gas bubbles, particularly for those consisting of dispersed gas bubbles (i.e., no considerable bubble stream). The gas plumes observed here should represent high gas flux and the rapid discharge of free gases. The maps of the backscattering intensity show that some bright patches scatter around the seabed (Figures 5, 6). The enhancement of the intensity should be attributed to the presence of authigenic carbonates accompanied by the free methane gases in the sub-seafloor sediments (Fonseca, 2001;Naudts et al., 2008). The cores of the authigenic carbonates associated with methane seepage have been got near Db4 (D5 in Figure 5), and anomalous backscattering intensity has also been observed there. The bicarbonate required for the formation of these carbonates is from the anaerobic oxidation of methane (AOM, Dale et al., 2008), which suggests high-flux methane had passed through the areas of these carbonates for years. Therefore, the bright patches of the backscattering intensity indicate the places of gas escaping that were active in the recent past.
The locations of the positive reliefs have a good match with those of the patches of the medium-to high-backscattering intensity (Figures 5, 6). Some hypotheses have been proposed how the sub-conical cones formed. They could result from the vertical displacement of sediments during the expulsion of pore fluids, the buoyancy of vertically continuous gas accumulation (Koch et al., 2015), or volume expansion during the formation of gas hydrate pingoes (Serié et al., 2012). These mechanisms point toward the occurrence of free gases bypassing the shallow underconsolidated sediments, though which one works here is not clear. The Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 557634 positive reliefs and the authigenic carbonates can both be taken as the geological expressions of the intermittent methane-rich gas emissions. These expressions can sometimes coalesce at the level of the seabed to form a larger one (e.g., ridge-like positive relief, Figures 5, 6).

Gas Emission in Transtensile Regime
The close spatial relationship between the faults and the acoustic features suggests the role of the faulting in determining the potential locations of gas emission. Consistent with the interpretations in the previous studies (e.g., Kimura, 1985;Gungor et al., 2012), the NE-SW oriented faults here are interpreted to form during the episodic back-arc rifting of the Okinawa Trough. It has to be stated that their trends are slightly different from those observed at the seabed in the deep-water setting ∼45 km to the east of the study area ( Figure 2). The difference of the fault trend between the deep-water setting and the slope here could be explained by the oblique rifting in the northern and the middle part of the Okinawa Trough (Fabbri et al., 2004). Its physical modeling results confirm this difference between the rift axis zone and the slope (Autin et al., 2013). It is very likely that part of the free gases in the subsurface escaped along the vertically permeable sediments of the damage zone of the faults that formed by back-arc rifting. The resultant extension would affect gas emissions more in the area closer to the rifting center ( Figure 1; Li et al., 2021), in which the strike-slip component is few. Therefore, the faulting driven by back-arc extension provided the migration pathways and the preferential exits for the free gases, and thus determines the potential locations of gas emissions at the seabed. Apart from the tensile stress, the shear one also has a contribution to the gas emissions ( Figure 10). The start of the gas emission requires lasting building up of the pore pressure and this can be primed by the local compression in the transtensive settings. The excess pore pressures are expected to occur in the restraining bend ( Figure 10). This is supported by the numerical modeling results showing that the relatively high σ 1 within the bend (Figure 9). These modeled stresses are almost perpendicular to that of the segment between the fault tips (red dashed line in Figure 9) and cause local compression toward it. This compression has also been evidenced by sediment thickening close to the restraining bend in a sand-box model (McClay and Bonora, 2001). Overall, the shallow sediments near these reliefs are most likely to be compressed as the tectonic activity continued and the resultant increasing pore pressure primed the formation of the gas emissions ( Figure 10). How these gas emissions formed is different from that in the Lusi mud volcano in Indonesia. The gas emissions there are located very close to a planar strike-slip fault and its activity is interpreted to significantly reduce the critical fluid pressure required to induce sediment fluidization (Mazzini et al., 2009). It is concluded that the restraining bend, other than the strike-slip fault planes, could be the potential place for the gas emissions bypassing the shallow marine sediments.

CONCLUSION
Multiple gas emissions are distributed unevenly along the faults at the seabed in the western slope of the Mid-Okinawa Trough. Part of the gas bubbles formed the gas plumes after escaping from the seabed, but did not reach the shallow seawater. Some other gases were trapped immediately below the domed seabed and at the places where there are no local morphological features. As a result of the ongoing back-arc rifting, the faulting allowed the free gases in the subsurface to pass through the shallow sediments. The components of the shear stresses caused local compression within the bending zone of the fault and the resultant increased pore pressure primed the formation of the gas emissions.

AUTHOR CONTRIBUTIONS
AL offered the key idea, designed and ran the FEM modelling, and drafted the manuscript. FC, NW and QL managed the project of the geophysical survey, provided the available data and the context of the methane seeps, vents and emissions in the Okinawa Trough. GY, GD and DL did the work of data processing and visualization. YS and XW checked the validity of the FEM model and made comments about the manuscript text.