Detection of Gas Hydrates in Faults Using Azimuthal Seismic Velocity Analysis, Vestnesa Ridge, W‐Svalbard Margin

Joint analysis of electrical resistivity and seismic velocity data is primarily used to detect the presence of gas hydrate‐filled faults and fractures. In this study, we present a novel approach to infer the occurrence of structurally controlled gas hydrate accumulations using azimuthal seismic velocity analysis. We perform this analysis using ocean‐bottom seismic data at two sites on Vestnesa Ridge, W‐Svalbard Margin. Previous geophysical studies inferred the presence of gas hydrates at shallow depths (up to ~190–195 m below the seafloor) in marine sediments of Vestnesa Ridge. We analyze azimuthal P‐wave seismic velocities in relation with steeply dipping near‐surface faults to study structural controls on gas hydrate distribution. This unique analysis documents directional changes in seismic velocities along and across faults. P‐wave velocities are elevated and reduced by ~0.06–0.08 km/s in azimuths where the raypath plane lies along the fault plane in the gas hydrate stability zone (GHSZ) and below the base of the GHSZ, respectively. The resulting velocities can be explained with the presence of gas hydrate‐ and free gas‐filled faults above and below the base of the GHSZ, respectively. Moreover, the occurrence of elevated and reduced (>0.05 km/s) seismic velocities in groups of azimuths bounded by faults suggests compartmentalization of gas hydrates and free gas by fault planes. Results from gas hydrate saturation modeling suggest that these observed changes in seismic velocities with azimuth can be due to gas hydrate saturated faults of thickness greater than 20 cm and considerably smaller than 300 cm.


Introduction
Structural and stratigraphic features play a significant role in controlling dynamic fluid flow processes occurring in the subsurface (Bjørlykke, 2015;Spencer, 2012). Low-density fluids, like hydrocarbon gases, tend to move upward from deeper to shallower depths due to pressure differences unless some impermeable stratum or structural feature traps them. These fluid flow processes can result in seepage of hydrocarbon gases from the seafloor (King & MacLean, 1970) and can also lead to shallow gas accumulations (Dondurur et al., 2011;Judd & Hovland, 1992;Vadakkepuliyambatta, 2014). In petroliferous basins, the distribution of fluids depends on fluid migration pathways within a basin, and thus, these pathways are very important to study (Caine et al., 1996;Davies & Handshy, 2003;Ligtenberg, 2005). Different structural features, like the presence of faults, fractures, and structural traps, determine migration pathways and accumulation zones for fluids (Bjørlykke, 2015;Suess et al., 2013;Vadakkepuliyambatta et al., 2013).
Faults play an important role in fluid migration systems, as these can be sealing or can also act as conduits, depending on the material present within a fault plane and the regional stress regime to which they are subjected (e.g., Caine et al., 1996;Davies & Handshy, 2003;Ligtenberg, 2005;Sibson, 1994). Hydrocarbons form at different depths depending on the process leading to their formation. Hydrocarbons formed from thermal degradation of organic matter, or by abiotic processes, generally have deeper origins compared to hydrocarbons formed by microbial degradation of organic matter (Etiope & Lollar, 2013;Schoell, 1988). Hydrocarbon gases migrate from deeper to shallower depths through permeable pathways, and in some cases, these migrated gases along with in situ produced biogenic gases get locked in shallow marine sediments in the form of gas hydrates (Kvenvolden, 1993).
Gas hydrates form under low temperature and high-pressure conditions and have gases trapped inside water molecules in a crystalline structure (Sloan, 1998). Gas hydrates occur widely in continental margins, primarily (99%) as methane hydrates (Ruppel & Kessler, 2017). Gas hydrates are stable only up to a certain depth below the seafloor, as temperature is too high to be compensated by pressure after this depth. A bottomsimulating reflection (BSR) often occurs at the base of the gas hydrate stability zone (GHSZ) as free gas is often trapped below gas hydrate-saturated sediments with low permeability (Shipley et al., 1979).
Apart from pressure and temperature requirements, the input of hydrocarbon gas and the presence of enough pore water are prerequisites for the formation of gas hydrates. Migration pathways can control the availability of these fluids; hence, gas hydrate distribution is affected by the capacity of structural features to accumulate and transport gas (Ben-Avraham et al., 2002;Bünz et al., 2012;Hennes et al., 2004;Taylor et al., 2000).
Pressure cores and well logs are well suited to investigate small-scale heterogeneities in the distribution of gas hydrates (Cook, 2010). These methods however, are not suitable for regional studies as they are restricted to individual stations. Seismic velocity analysis provides a means to study the distribution of gas hydrates (Chand et al., 2004;Kumar et al., 2006;Madrussani et al., 2010). Seismic velocities estimated using oceanbottom seismic (OBS) data have been successfully used in different geological settings to estimate gas hydrate saturations in the GHSZ (Bünz et al., 2005;Kumar et al., 2007;Satyavani et al., 2013;Song et al., 2018;Westbrook et al., 2008). Gas hydrate bearing sediments have a higher bulk modulus and thus exhibit higher seismic velocities (Chand et al., 2004;Ecker et al., 1998;Gei & Carcione, 2003;Lee & Collett, 2001). Changes in seismic velocities can potentially indicate pore-fill variations, that is, the presence of free gas, gas hydrates, or perhaps the presence of authigenic carbonates where fluid focusing occurs (Kumar et al., 2006;Pecher et al., 2003;Toksöz et al., 1976). The identification of small-scale (~10-200 m) lateral changes in seismic velocities due to variations in distribution of gas hydrates in pore spaces and faults/fractures is challenged by limitations in the seismic resolution (Kumar et al., 2006). Thus, seismic velocity analysis is rarely used to investigate small-scale lateral changes in gas hydrate and free gas saturations (Jaiswal et al., 2012b;Satyavani et al., 2013). In the present study, we document results of a high-resolution OBS experiment where we correlate azimuthal P-wave velocity variations with fine-scale structural maps from 3-D seismic data. We aim to investigate structural controls, in proximity to fault structures, on the gas hydrate and free gas distribution, and on the focused fluid flow in Vestnesa Ridge, west-Svalbard margin.

Study Area
Vestnesa Ridge is a large sediment drift located at~79°N on the western Svalbard continental margin (Eiken & Hinz, 1993;Howe et al., 2008;Hustoft et al., 2009; Figure 1). It bends from southeast-northwest direction to east-west direction toward the north (Figure 1b). It is bounded by the Knipovich Ridge in the south, Molloy Ridge in the northwest and Molloy Transform Fault in the southwest (Figure 1a). Rifting at these mid-ocean ridges (Knipovich Ridge and Molloy Ridge) and shear motion at Molloy Transform Fault dictate regional stresses and faulting patterns in the region (Plaza-Faverola & Keiding, 2019). Sedimentation on the western Svalbard margin occurred under the influence of bottom-water contourite currents and mainly consists of turbiditic, glaciomarine, and hemipelagic sediments (Eiken & Hinz, 1993;Ottesen et al., 2005;Stein et al., 2005).
The presence of a gas hydrate system in Vestnesa Ridge is well documented (Bünz et al., 2012;Goswami et al., 2015;Hustoft et al., 2009;Petersen et al., 2010;Singhroha et al., 2019;Smith et al., 2014;Vogt et al., 1994). Seafloor pockmarks and gas chimneys along the ridge are linked to faults and fractures (Plaza-Faverola et al., 2015;Singhroha et al., 2016). Considering the suggested presence of thermogenic gas (Panieri et al., 2017;Plaza-Faverola et al., 2017;Smith et al., 2014) in the study area, the study of fault/fracture systems is even more important as thermogenic gases have deeper origins (often below the GHSZ) and structural features affect the migration pathways of thermogenic gases. Faults at deeper depths in this region can be potential fluid migration pathways for deep sourced warm fluids (Knies et al., 2018;Dumke et al., 2016;Waghorn et al., 2018), whereas faults at the shallower depths within the GHSZ can be migration pathways or can also act as seals, as they can potentially be plugged with gas hydrates (Goswami et al., 2017;Madrussani et al., 2010).
Tectonic stresses also play an important role in deciding whether a fault is sealing (blocks the passage of fluids) or nonsealing (allows the passage of fluids). Stress field variations and the kinematics of faults through geological time have been suggested to exert a control on the temporal and spatial distribution of seepage along Vestnesa Ridge (Plaza-Faverola et al., 2015;Plaza-Faverola & Keiding, 2019). The eastern segment of Vestnesa Ridge has several pockmarks, through which methane seeps, compared to the western segment, where pockmarks are inactive (Bünz et al., 2012). This is potentially due to the alignment of tectonic stresses, which govern the opening and sealing of faults in the eastern and western segment of Vestnesa Ridge, respectively (Plaza-Faverola et al., 2015;Plaza-Faverola & Keiding, 2019). P-and S-wave velocity analysis (Singhroha et al., 2019) and Q analysis (Singhroha et al., 2016) suggest that faults control the distribution of gas hydrate and free gas in the eastern segment of Vestnesa Ridge. Hence, previous studies suggest that the structural settings from regional (>30 km) to local scale (1-2 km) dictate the distribution of fluids within the gas hydrate system.

Structural Evolution of Faults in Vestnesa Ridge
Vestnesa Ridge has several large-and small-scale fault features. The eastern segment of Vestnesa Ridge has more faults in the subsurface compared to the western segment of Vestnesa Ridge (Plaza-Faverola et al., 2015). Seismic variance attribute maps indicate the presence of fine-scale (a few meters resolution) faults and fractures associated with gas chimneys and seafloor pockmarks (Plaza-Faverola et al., 2015). These maps give a good overview of the evolution of subsurface faults over time (Figures 1e-1g). We pick three horizons in the eastern segment of Vestnesa Ridge ( Figure 1d) and estimate the seismic variance attribute along these three horizons to study the evolution of faults and fractures.
The northeastern flank of the eastern segment of Vestnesa Ridge has large-scale (~6-10 km) extensive faults (for example, Fault3 in Figure 1g) typically associated with sediment slides. These faults are continuous even at deeper depths (Figures 1e and 1f). Faults in the southwestern flank of the eastern segment of Vestnesa Ridge appear relatively laterally smaller, particularly at shallow depths (for example, Fault1 and Fault2 in Figures 1e-1g). At shallower depths, these faults are spatially limited within close proximity of fluid flow features (Figure 1g), suggesting that high-fluid pressure could play a role in stress build-up and faulting patterns. Such changes might imply that high-fluid pressure triggers faulting in this region. At deeper depths, faults can act as potential pathways for fluid migration and fluids tend to migrate toward gas chimneys at shallower depths, thus making fault systems spatially limited at shallower depths. Such patterns are not observed in faults (for example, Fault3) on the northeastern flank of the ridge.

Seismic Data Acquisition
In the present study, we aim to analyze the role of faults on the distribution of fluids in the eastern segment of the Vestnesa Ridge gas hydrate system on a very small-scale (<200 m). We selected OBS locations near the Vestnesa Ridge crest and close (<30 m) to fault planes of Fault1 and Fault2 in survey1 and survey2, respectively (Figures 1c and 1e). Seismic velocity analysis (Singhroha et al., 2019) and seismic Q analysis (Singhroha et al., 2016) predicted differences in gas hydrate and free gas saturations across Fault2. Fault1 compared to Fault2 has relatively less lateral extent, spatially limited to only one side of the OBS location (Figures 1c and 1e). We performed seismic velocity analysis to study the detailed variation of seismic velocities with azimuth at these two OBS sites (Figure 1).
We acquired OBS data in a circular geometry, where we run the seismic source in concentric circular paths around two OBS locations (Figure 1c). We shot in three and four concentric circles in survey1 (in 2015) and survey2 (in 2016), respectively. Seismic energy was generated by a mini generator-injector (GI) airgun (15/15 in 3 ; Sercel) and a GI airgun (45/45 in 3 ; Sercel) in survey1 and survey2, respectively. The mini GI airgun provides higher resolutions compared to the GI airgun as the mini GI airgun generates a source signal that has broader amplitude spectrum (10-300 Hz) and contains significantly higher energy in high frequencies (peak frequency~150-180 Hz). High frequencies produced by mini GI and GI airguns make reflections sharp, allowing reflections to be picked with a low pick uncertainty (<1 ms). Seismic data were sampled at every 0.25 ms. Data processing included standard band pass filtering to improve the signal-to-noise ratio and to facilitate picking of seismic horizons.
This circular acquisition approach has been used earlier for tomographic seismic velocity and seismic anisotropic analysis by qualitatively looking at changes in seismic amplitudes and arrival times in different directions (e.g., Exley et al., 2010;Haacke & Westbrook, 2006;Plaza-Faverola et al., 2010;Satyavani et al., 2013). The directional seismic velocity information is lost in seismic velocity analysis from the tomographic approach, as seismic velocity is assumed the same for a ray within a grid, regardless of the direction it enters the grid. We use this acquisition approach further to estimate azimuthal seismic velocities, by picking travel times in different azimuths and then inverting these travel times.

Detection of Gas Hydrates in Faults and Fractures
The presence of gas hydrates in faults and fractures is primarily detected using resistivity methods, as resistivity estimates are very sensitive to the occurrence of gas hydrates in faults and fractures (Cook, 2010). Seismic velocity estimates, assuming an isotropic and homogeneous medium, are comparatively far less sensitive to the occurrence of gas hydrates in faults and fractures (Cook, 2010;Cook & Waite, 2018;Ghosh et al., 2010;Lee & Collett, 2012). Differences observed in gas hydrate saturation estimates, from resistivity and seismic velocity analysis, potentially indicate the occurrence of gas hydrate-filled faults and fractures (Cook, 2010;Singhroha et al., 2019).
Often faults and fractures in a region are aligned in a direction, dictated by local and regional stresses (Bayly, 1992). However, faults and fractures generated by fluid-overpressure may not follow the regional stress field (Bayly, 1992;Cox, 1995;Sibson & Scott, 1998). Vestnesa Ridge may potentially have faults and fractures generated by both mechanisms (Plaza-Faverola et al., 2015). Some of the faults in the study area are spatially limited and often have smaller fault/fracture networks popping out of a prominent fault (Figures 1e-1g). With the circular seismic acquisition experiment, we can illuminate these structures, which allow us to study the properties of the material filling in the deformation planes ( Figure 2). For instance, fast seismic velocities for seismic rays whose raypath plane lie along a fault plane, hint that the secondary porosity has been infilled with a fast seismic velocity material ( Figure 2). The schematic diagram shows raypaths arriving from different azimuths. The violet color shows rays with raypaths passing through fault and fractures in the subsurface. The green color shows rays arriving from the back side of the fault plane, whereas the blue color shows rays passing arriving from the front side of the fault. The green and blue color rays have raypaths that do not pass through fault or fractures. Haacke and Westbrook (2006) studied converted shear waves from near-surface sediments on the West-Svalbard Margin and found changes in the particle motion and the polarization direction with depth suggesting the presence of seismic anisotropy in the shallow sub-surface (approximately 200-300 m below the seafloor). They associated these changes with the presence of gas hydrates in aligned cracks and suggested the study of detailed azimuthal anisotropy in the area.
Properties of a medium change with the direction of wave propagation in an anisotropic medium (Crampin, 1981). This occurs due to granular (or crystal) anisotropy, leading to microscopic directional variations in seismic wave propagation properties in some minerals (for example, quartz and clay) or due to preferential alignment of faults/fractures, minerals, and rocks in a given direction (Musgrave, 1970;Nelson, 1985). Changes in mineral/rock composition or distribution in a medium create heterogeneity in a medium. Every heterogeneity also generates anisotropy in a very broad sense or definition of anisotropy (Winterstein, 1990). However, in a strict sense, we consider medium anisotropic only if the heterogeneity generating anisotropy is smaller than the wavelength of seismic waves (Winterstein, 1990). Hence, changes in seismic properties due to the presence of faults and fractures (with fault thickness considerably smaller than seismic wavelengths) are considered anisotropic changes in a medium. Vertical faults and fractures create azimuthally anisotropic media and affect the propagation of seismic waves (Zheng, 2006). However, the presence of spatially limited faults and fractures may limit the application of anisotropic principals to this study as anisotropic studies assume the presence of infinitely laterally extended faults in the horizontal direction (Winterstein, 1990). In addition, seismic rays from different directions should pass through a same geologic body to analyze seismic anisotropy in its strict sense. The acquisition geometry in this study does not allow us to illuminate the same geologic body from different directions due to fixed receiver locations.
Faults and fractures are often filled with minerals and fluids that are different from the background matrix creating a notable difference in their bulk effective properties. We expect changes in seismic velocities with azimuth due to changes in the degree of illumination of faults and fractures by seismic raypaths with azimuth. By changes in illumination, we mean variations in raypath lengths travelling through potentially high-seismic velocity (for example, gas hydrate) or low-seismic velocity (for example, free gas) material filled in faults and fractures ( Figure 2). Hence, we attempt to study seismic velocity properties of material filled in faults and fractures by acquiring seismic data in a circular geometry around OBS sites.
In a circular seismic acquisition geometry, it is easy to pick reflection arrivals at different offsets in different azimuths, therefore an appropriate approach for the study of azimuthal variations in seismic velocities ( Figure 3). Travel times from two different offsets corresponding to a reflection from a given flat layer boundary are ideally sufficient to get two unknown parameters in a layer, that is, the thickness and its seismic velocity. Thus, travel times corresponding to three and four different offsets in different azimuths in survey1 and survey2, respectively, are theoretically enough to constrain the subsurface seismic velocity model. A simple flat-layered model with five units has been used to test the validity of this approach (Figure 4a). We use the layer stripping approach to constrain the seismic velocity and the thickness of the fourth layer (Figure 4a) using reflection arrival times at different offsets. Different subsurface seismic velocity models can give the same reflection arrival time at a given offset. Each curve in Figures 4b and 4c shows different possible combinations of seismic velocity and thickness for the fourth layer ( Figure 4a) that will have the same reflection arrival time at the given offset. When thickness-seismic velocity pairs must accommodate travel times for different offsets, the number of possible solutions is significantly reduced, which increases the uniqueness of the seismic velocity model. Figures 4b and 4c show the convergence of curves from different offsets at a single point, that is, at the true seismic velocity and thickness of the layer. Travel times from different offsets, especially far offsets, are needed to increase the confidence in the derived seismic velocity model. For arrival times corresponding to nearby offsets, possible seismic velocity-thickness curves are very close, and seismic velocity and thickness values obtained by the intersection of these curves will produce results that will have a very low confidence. Therefore, if few accurate arrival times can be picked at near, medium, and far offsets; reliable seismic velocity models can be derived as curves for different possible seismic velocity, and thickness parameters will intersect each other at very high angles that will produce a seismic velocity model with high confidence (Figures 4b and 4c). We also test the sensitivity of the seismic velocity estimates from this approach due to the potential uncertainty in the picked reflection arrival times. Results show that an error of 1 ms in the arrival times corresponds to a maximum possible seismic velocity error of 0.035 km/s in the final velocity model (Figure 4c).

Seismic Velocity Modeling
In order to derive an azimuthal seismic velocity model, seven prominent reflection arrivals were picked at both investigated sites in a seismic section acquired using circular acquisition geometry (azimuth and offset vary with every trace in the seismic section as shown in Figure 3a). These picked reflection arrival times were sorted by the azimuth and the offset between shot points and the OBS instrument. A seismic velocity model in a given azimuth was calculated using travel times of different reflection arrivals along this azimuth (Figures 3b and 3c).
We use high-resolution 3-D P-Cable seismic data (Plaza-Faverola et al., 2015; Figure 1c) and OBS data (Singhroha et al., 2019) acquired in this area to get a starting model for azimuthal seismic velocity analysis. Seismic velocity models estimated using OBS data (Singhroha et al., 2019) are used to convert the 3-D P-Cable seismic data in the depth domain. Two-dimensional seismic velocity models for different azimuths are then extracted from the depth-converted 3-D seismic data. These 2-D models provide very good estimates of dips in dipping layers and constrain the geometry of layers in the region. With these 2-D models as initial seismic velocity models, we use Rayinvr program based on Zelt and Smith (1992) approach to invert the picked travel times in a layer stripping technique to derive the final seismic velocity model. The travel time inversion results for picked travel times along 10 and 190°azimuths are shown in Figures 3b and 3c. Similar analysis is done in 36 different directions each containing 10°radial azimuths to get an overview of variation of velocities with azimuth around the OBS stations (Figures 3, 5, and 6). The root-mean-square misfit between the picked arrival times and travel times estimated using Rayinvr in the final seismic velocity model is less than 1 ms, therefore indicating a great degree of confidence in the derived seismic velocity models (Figure 3). A low pick uncertainty (<1 ms) and a low root-mean-square misfit (<1 ms) in the best fit seismic velocity models ensure a low seismic velocity uncertainty (~0.035 km/s as shown in Figure 4c) in azimuthal seismic velocity estimates. Azimuthal seismic velocities are estimated for layers L1, L2, L3, L4, L5, and L6 in survey1 (Figures 3a and 5) and for layers LL1, LL2, LL3, LL4, LL5, and LL6 in survey2 ( Figure 6). We attempted to pick the same layers in both sites. However, the difference in data resolution (due to differences in seismic source as discussed in section 4) and the change in strength and discontinuity of the picked phases at two sites prevented us from picking same phases. A distinct prominent BSR can be seen in both the sites, and it lies in survey1 between L5 and L6 and in survey2 between LL5 and LL6.

Results From Azimuthal Seismic Velocity Analysis
Azimuthal seismic velocities for each individual layer estimated using travel time inversion are plotted using radial pie charts, with different directions representative of different azimuths and the radius of the circle representative of the spread of ray path (Figures 5 and 6). In layer L1 of survey1 (layers shown in Figure 3

10.1029/2019JB017949
Journal of Geophysical Research: Solid Earth a), the variation in seismic velocity is small (1.506-1.527 km/s) and falls within the uncertainty (0.035 km/s) limits ( Figure 5). In layer L2, a seismic velocity change is observed from around 1.57 km/s toward the southwest to around 1.63 km/s within an azimuthal fan of around 60°in the east-northeast region ( Figure 5). Small elevations in seismic velocities (0.015-0.020 km/s) in 10°azimuthal fans compared to the average velocity (around 1.61 km/s) are also observed in azimuths corresponding to 300 and 320°( Figure 5). In layers L3 and L4, seismic velocities in one-half of the azimuths toward the northeast (~315-135°) are~0.03-0.04 km/s higher compared to the average velocity (~1.62-1.64 km/s) in the other half ( Figure 5). In layer L5, high seismic velocities (1.84-1.86 km/s) occur in the north-northeast direction. A patch (20°radial azimuthal fan) of a particularly high seismic velocity (>1.87 km/s) is observed in the northwest direction. Seismic velocities are lower (~1.79-1.80 km/s) toward the southeast ( Figure 5). In layer L6 (layer below the BSR), seismic velocities are considerably lower than in the layers above ( Figure 5). The average seismic velocity in this layer is 1.343 km/s, and it is relatively uniform in different azimuths, except in few orientations (10°radial azimuthal fans along 100°, 110°, 190°, and 310°azimuths) where seismic velocities are slightly lower ( Figure 5).
Radial pie charts are also used to plot results from survey2 ( Figure 6). In layers LL1 and LL2, high seismic velocities (1.54-1.56 km/s) are observed in the northeast direction compared to the average velocity (~1.50-1.52 km/s). Layer LL3 shows an average seismic velocity of 1.68 km/s, with two azimuthal fans of slightly elevated seismic velocities at 30°(1.74 km/s) and 120°(1.714 km/s). In layer LL4, a patch of relatively higher seismic velocity (1.73-1.745 km/s) compared to an average velocity of 1.71 km/s in this layer lies in the south-southeast direction. In layers above (LL5) and below (LL6) the BSR, very strong variations in seismic velocities can be observed in a number of azimuths ( Figure 6b). In layer LL5, seismic velocities on the southwestern side (average velocity around 1.76 km/s) are higher than seismic velocities on the northeastern side (average velocity of 1.73 km/s). Higher seismic velocities are observed along the radial fan of around 40°i n the southeast direction with seismic velocities reaching up to 1.823 km/s. The average seismic velocity in layer LL6 (1.41 km/s) is considerably lower. In LL6, particularly low seismic velocities (1.36-1.38 km/s) are observed in a radial fan of 40°in the southeast direction and along 320°azimuth. In the azimuths corresponding to northwest directions, relatively higher seismic velocities (~1.45-1.50 km/s) are observed (Figure 6b).

Correlation Between Seismic Velocity Anomalies and Structural Features
We analyze our azimuthal seismic velocity results along with seismic variance maps that show structural trends in the study area (Figures 1e-1g, 5, 6c, 7, and 8). For some azimuths, the raypath plane lies along the direction of the fault plane and significant raypath lengths traverse through that fault plane in the subsurface (as shown in Figures 2 and 9c). The overlay of seismic velocities with seismic variance maps allows us to correlate seismic velocity anomalies with faults and fractures in the region (Figures 5-8). Horizon H30 (Figure 1d) that lies around 30 m below the seafloor is used to analyze anomalies in the shallow layers ( Figure 7). Two deeper horizons, that is, H50 and H80 (Figure 9) corresponding to approximately 0.3 and 1.5 Ma, respectively (Plaza-Faverola et al., 2015), are used to interpret anomalies for layers that lie in the middle (around 70-m depth below the seafloor; Figure 8) and near the base of the GHSZ (Figures 5 and  6), respectively.
In survey1 and survey2, a very strong relationship is observed between the distribution of seismic velocity anomalies and faults in the subsurface. At shallow depth (<40 m below the seafloor; Figure 7), the link between faults and anomalies in azimuthal seismic velocity is evident (Figure 7). In survey1, there are hardly any faults as indicated by a smooth seismic variance map (Figure 7). Such a smoothness correlates with a largely invariable (1.50-1.53 km/s) azimuthal seismic velocity (Figure 7). The opposite is true for survey2 where seismic velocities vary (1.49-1.56 km/s) and elevated velocities (1.55-1.56 km/s) in the northeast direction seem to be related to the fault (marked in Figure 7) in the northeast direction. Seismic velocity anomalies are confined by this fault and another fault in the area (Figure 7). There is also a small increase (~0.02 km/s) in the seismic velocity in south-southeast direction (i.e., in the direction of the marked fault in Figure 7).
The correlation between azimuthal seismic velocities and structural features is even clearer for the seismic variance map from H50 horizon (Figure 8). In survey2, seismic velocities vary (1.64-1.74 km/s) and elevated (1.74 km/s) seismic velocity in the northeast direction seems to be related to the fault (Figure 8 Figure 8). In addition, there are high (>1.62-1.63 km/s) seismic velocities in the northeastern half toward the ridge crest, except in the shadow zone (i.e., azimuths whose raypaths lie northeast to the fault), where low seismic velocities occur (<1.61 km/s; Figure 8).
Changes in seismic velocities along and across faults as observed at horizon H30 and H50 depths give an indication about the effect of faults on seismic velocities (Figures 7 and 8). However, these seismic velocity changes (~0.015-0.05 km/s) are often within or around the uncertainty limit (0.035 km/s). The relationship between faults and azimuthal seismic velocity is most evident near the BSR where changes in seismic velocities due to faults are well above the uncertainty limit (Figures 5 and 6). In layer L5 (the layer above the BSR), a spike (~0.06-0.08 km/s) in the seismic velocity in the northwest direction clearly seems to be related to the presence of a fault ( Figure 5). There is also a decrease (~0.04 km/s) in the seismic velocity in the shadow zone (i.e., azimuths whose raypaths lie northeast to the fault; marked as S in the azimuthal seismic velocity plot for L5 in Figure 5). In layer L6, a decrease (~0.04 km/s) in seismic velocity occurs along the fault ( Figure 5). In layer LL5 (layer above the BSR) in survey2, where the OBS is placed near Fault2, seismic velocity spikes (~0.06-0.08 km/s) are observed in between faults shaping an X (Figures 6b and 6c). The orientation and alignment of faults also affect the regional distribution of seismic velocities (Figures 6b and 6c). Seismic velocities in the northeast direction are generally lower than velocities in the southwest direction. In layer LL6 (the layer below the BSR), a decrease (~0.06-0.08 km/s) in the seismic velocity occurs along

Occurrence of Gas Hydrate-Filled Faults
Directional illumination of subsurface by seismic rays followed by the estimation of azimuthal seismic velocities can be a potential technique to detect the material filled in faults and fractures. The striking correlation of seismic velocity variations and apparent anomalies with the orientation and position of subsurface structures suggests that the approach has been implemented successfully. Small azimuthal seismic velocity variations (<0.02 km/s), which are well below uncertainty limits (<0.035 km/s), especially in deeper layers, can be due to small errors (~1 ms) in travel-time picks and differences in raypaths in shallower layers ( Figure 9). However, large changes in seismic velocities (>0.03-0.04 km/s) observed along azimuths, where significant length of raypaths fall in fault planes (Figure 9c), may be due to material filled within faults in the region. The correlation of seismic velocity anomalies with faults is stronger with depth (especially above and below the BSR). This correlation is even more prominent in survey2 that has more laterally extended and better defined fault planes passing through the OBS location compared to survey1 that has a relatively small-scale fault system (Figures 1c-1g, 5, and 6). In addition, Fault1 in survey1 extends only on one side of the OBS location and thus gets illumination by seismic raypaths from a limited group of azimuths in the northwest direction. (Figures 1c, 1e, and 5). However, Fault2 in survey2 laterally extends across the OBS location and thus gets bimodal illumination from a wider group of azimuths in the eastward and 180°o pposite westward direction (Figures 1c, 1e, and 6b).
The observed azimuthal seismic velocity contrasts therefore indicate whether secondary permeability created at damage zones nearby faults has favored the accumulation of any fluid. Variations in fluid content and concentration has a big impact on the bulk properties of the medium (Batzle & Wang, 1992; Brown & Figure 9. Schematic of raypaths for L5 in Inline and crossline directions at the survey site 1 (9a and 9b). (c) Schematic of raypaths for LL5 arriving from three different azimuths (with four rays for each azimuth) along with the seismic variance map at the survey site 2. The green color shows raypaths of downgoing rays, whereas the dotted red color shows raypaths of upgoing rays (reflected from the boundary between layers LL5 and LL6). The angle between the fault plane and raypath plane is very low for four rays (four rays on the left side), which corresponds to one azimuth. Length of raypaths passing through the fault plane will be higher for these set of rays compared to rays from other azimuth whose raypath plane intersect at high angle with the fault plane. Faults are near vertical in this area. Korringa, 1975;Gassmann, 1951;Han, 1992;Mavko et al., 1995Mavko et al., , 1998Nolen-Hoeksema, 2000;Sengupta & Mavko, 1999). Secondary porosity, in the form of faults and fractures, plays an important role in storing and transporting fluids and different fluids have significant impact on the seismic velocity characteristic of a medium (Boadu & Long, 1996;Brown & Scholz, 1986;Hudson, 1981). Singhroha et al. (2019) and Goswami et al. (2015) document a gas hydrate system in this area where different components of the gas hydrate system, that is, methane hydrate and free gas, strongly affect the bulk seismic velocity. The presence of gas hydrates increases the P-wave velocity, and the presence of free gas in sediments decreases the P-wave velocity (Chand et al., 2004;Ecker et al., 1998;Gei & Carcione, 2003;Lee & Collett, 2001;Song et al., 2018). High seismic velocity anomalies observed in the GHSZ above the BSR in some azimuths (Figures 5-8) can be due to the presence of gas hydrates. Other potential high seismic velocity material like authigenic carbonates can be ruled out as they do not form at this depth and there is no indication of paleo seepage activity at the faults. Low seismic velocity anomalies observed in the layer below the BSR (Figures 5 and 6) is a classic indicator of free gas accumulations beneath gas hydrate bearing sediments (Goswami et al., 2015;Singhroha et al., 2019).
In survey1 and survey2, a seismic velocity increase (up to~0.06-0.08 km/s in layers L5 and LL5) in layers in the GHSZ with azimuths corresponding to fault orientations compared to other azimuths ( Figures 5-8) suggests a preferential accumulation of hydrates within faults. A seismic velocity decrease (~0.04 km/s in layer L6 and~0.06-0.08 km/s in layer LL6) within the free gas zone along the fault direction compared to other directions ( Figure 5) further proves the role of faults in directing gas migration toward the GHSZ (Bünz et al., 2012;Plaza-Faverola et al., 2015;Singhroha et al., 2016 ;Singhroha et al., 2019). In addition, near the BSR depth in survey2, increase (~0.06-0.08 km/s) and decrease (~0.06-0.08 km/s) of seismic velocities occur between fault planes above and below the BSR, respectively (Figures 6b and 61c). This potentially indicates compartmentalization of gas hydrates and free gas by fault planes above and below the BSR, respectively (Figures 6b and 6c).
The evidence of occurrence of gas hydrates in faults is stronger at deeper depths. In shallow depths, small changes (~0.015-0.05 km/s) in seismic velocities are observed along fault orientations. This difference can be because of two reasons. Firstly, the subsurface is relatively less faulted near the seafloor. Another reason can be the fact that there is lower gas hydrate saturation at shallow depths and higher gas hydrate saturation near the base of the GHSZ (Singhroha et al., 2019). Both these factors reduce the likelihood of the presence of gas hydrate-filled faults at shallow depths. In addition, these changes are close to or below the uncertainty limit (0.035 km/s), which decreases our confidence in associating seismic velocity changes at shallow depths with the presence of gas hydrate in faults.

Effect of Faults on Gas Hydrate and Free Gas Distribution
In survey1, higher seismic velocities obtained in one half of azimuths toward the northeast (~315°-135°) in layers in the GHSZ (L1 to L5; Figure 5) are consistent with higher concentrations of hydrates toward the ridge crest as reported based on multi-channel seismic, 2-D OBS surveys and 3-D seismic attenuation studies (Hustoft et al., 2009;Singhroha et al., 2016;Singhroha et al., 2019). However, elevated (~0.03-0.04 km/s) seismic velocities in layer L5 toward the downslope (southwest) side of Fault1 compared to seismic velocities in the shadow zone (marked as S in Figure 5), that is, toward the upslope (northeast) side of Fault1 show the effect of faults on the distribution of gas hydrates in the region ( Figure 5). Raypaths in these two cases lie close (<0.04 km) to Fault1 but on opposite (northeast and southwest) sides to this fault, and raypaths do not intersect this fault (Figures 5 and 9). Hence, changes in seismic velocities in azimuths across Fault1 potentially indicate changes in gas hydrate saturation across Fault1. A relative increase in seismic velocities in the northeastern half observed in layer L5 in survey1 ( Figure 5) is not observed in layer LL5 in survey2 (Figures 6b and 6c). On the contrary, on average, seismic velocities in layer LL5, along azimuths in the southwestern half of Fault2, are slightly (~0.02-0.03 km/s) higher than seismic velocities along azimuths in the northeastern half of Fault2 (Figures 6b and 6c). This can be due to extensive shadow zone (in upslope direction) created by laterally extended Fault2. This shadow zone was limited to the selected group of azimuths (marked as S in Figure 5) in survey1, as Fault1 in this area only extends to one side of the OBS location. The average decrease in seismic velocities above the BSR in the shadow zone indicates decreased gas hydrate saturations in shadow zones. This potentially shows the control of faults on the upslope migration of fluids. Gases and other fluids trapped below impermeable gas hydrate saturated sediments in a layer below the base 10.1029/2019JB017949

Journal of Geophysical Research: Solid Earth
of the GHSZ move upslope (within the same layer below the GHSZ) toward the ridge crest from deeper (>2 km) depths (Bünz et al., 2012;Hustoft et al., 2009;Singhroha et al., 2016). This upslope migration of fluid and gases can be potentially obstructed by Fault1 and Fault2, thus leading to decreased gas hydrate saturations in the shadow zones due to the lack of availability of gases to form hydrates. Earlier estimates also show decreased gas hydrate saturations in the northeast direction from Fault2 (Singhroha et al., 2019). This theory is further concurred by higher seismic velocities in the free gas zone (layer LL6 in survey2) potentially suggesting lower free gas saturations in the northeastern half (i.e., in the upslope direction to the fault; Figures 6b and 6c). Assuming that this theory is plausible, it can be inferred that faults are sealed and hindering the upslope upward gas migration toward the ridge crest. Alternatively, faults are rather dilated and favoring gas leakage through these faults, thus limiting the presence of free gas on the other side of the fault.
In the GHSZ, faults can act as seals either due to stress or due to the presence of impermeable gas hydrates. The likelihood of faults being sealed in this segment of Vestnesa Ridge is relatively low as methane gas seeps through several pockmarks in this area (Plaza-Faverola et al., 2015;Plaza-Faverola & Keiding, 2019). If a fault is sealing due to the occurrence of gas hydrates in faults, we expect a relatively uniform free gas saturation below the BSR across the fault, as gas hydrate-filled sealed faults cannot affect the free gas distribution below. However, we find higher free gas saturations on the downslope side to the fault (Figures 5 and 6).
Results from seismic velocity analysis in layers L5, LL5, L6, and LL6 suggest a relatively low gas hydrate (in the GHSZ) and free gas saturation (in free gas zone) in a direction upslope (i.e., northeast) with respect to northwest-southeast oriented faults (Figures 1c-1g, 5, and 6). This is further corroborated by the fact that seismic velocity and seismic Q analysis show differences in free gas saturations across the fault (Singhroha et al., 2016;Singhroha et al., 2019). In addition, we find highest free gas saturations in the faults in the free gas zone (Figures 5 and 6). These observations hint toward a system in which gas migrates upward through fault networks from the base of the GHSZ. These fault networks spatially shrink at shallower depths and terminate within a gas chimney, which actively seeps methane (as described in section 3).
Continuous migration of gas through faults can be due to high fluid pressure from below. Gas migrating upward under the base of the GHSZ seems to be compartmentalized by the presence of faults in the GHSZ (Figure 6c). The fluid pressure southwest of the faults can be very high compared to northeast of the fault due to difference in free gas saturation. This fluid pressure can create fractures and faults southwest of the fault, and gas hydrate accumulation would favorably occur in these fractures, which are embedded in fine-grained hemiplegic sediments. These sealed faults can lead to accumulation of free gas that may result in an increase in overpressure from the gas below and a reactivation of sealed faults in the GHSZ (Flemings et al., 2003;Hornbach et al., 2004;Kleinberg, 2005;Liu & Flemings, 2007). These faults are spatially limited and connected with gas chimneys at shallow depths and thus can release gas flowing through these reactivated faults. The overpressured gas flowing through faults can generate fractures that may lead to further gas hydrate deposits in these fractures (Kumar et al., 2006;Yan et al., 2017).

Characterization of Gas Hydrate-Filled Faults
The presence of gas hydrate in faults increases the P-wave velocity compared to water saturated faults where the P-wave velocity slightly decreases as water is a weak inclusion (Hudson, 1981). Faults of finite (>10 cm) thickness (opposed to assumed paper-thin fault planes in the area) filled with gas hydrates can generate the observed azimuthal seismic velocity anomalies in the GHSZ. There are several established relations between hydrate content, morphology, and seismic velocities (Chand et al., 2004;Ghosh et al., 2010;Helgerud et al., 1999;Jakobsen et al., 2000;Lee et al., 1996). However, considering the complexity of raypaths through faults in this study, we have used simple time-average equations to estimate the effect of fully gas-hydrate saturated faults on seismic velocities ( Figure 10).
In this modeling approach, we have closely analyzed parameters such as the thickness of a fault plane and the offset between the fault plane and location of receiving station (i.e., OBS in this study). The offset between the fault plane and location of receiving station can also theoretically compensate for its undulating geometry of a fault plane and its tilt with respect to vertical. Assuming this offset to be 0.005 km (which approximately corresponds to 1.44°angle between the fault plane and the raypath plane), we estimate that 20-cm thickness of a gas-hydrate saturated fault can generate a seismic velocity anomaly of 0.072 km/s (Figures 10a-10c). Similarly, an anomaly of 0.072 km/s can also be roughly generated by fractures (for example, 40-cm fracture thickness and 50% fracture density and 80-cm fracture thickness and 25% fracture 10.1029/2019JB017949

Journal of Geophysical Research: Solid Earth
density; Hudson, 1981). Results from this modeling approach also suggest that there will be changes in seismic velocities at a detectable level (>0.02 km/s) only if the raypath plane and the fault plane lie in the same direction or if they intersect at a very small angle (<10-20°; Figures 10a-10c). We also did not find Figure 10. Schematic of the effect of gas-hydrate saturated faults on P-wave velocity. Plan (a) and side (b) views show a typical propagation of P-wave through a layer. P-wave passes through gas hydrate saturated fault as shown in different views. The effect of gas hydrate-saturated fault on P-wave velocity depends on the azimuth as length of raypath passing through gas hydrate-saturated faults vary with azimuth (c). Angle between the fault and the raypath plane (ф) decides the length of raypaths passing through a fault plane and thus different gas hydrate-saturated fault thicknesses at different ф can generate similar (0.07 km/s) seismic velocity anomaly (d).

10.1029/2019JB017949
Journal of Geophysical Research: Solid Earth any significant changes in seismic velocities due to a fault that have a fault plane aligned at a high angle (>20-30°) to raypath plane (for example, the fault shown in green color over azimuthal seismic velocity model for layer LL5 in Figure 6b). This is mainly due to changes in the angle between the raypath plane and the fault plane, which effectively changes the length of the raypath passing through a fault plane with azimuth (Figures 9c and 10a-10c).
Considering the observed seismic velocity spikes of around 0.06-0.08 km/s in azimuthal seismic velocity estimates, we argue that this can be due to a fully saturated gas hydrate fault of approximately 20-cm thickness or more in cases if fractured system coexits instead of fully developed gas hydrate-saturated fault system. These estimates assume the minimum angle between the fault plane and the raypath plane to be 1.44°. However, if the minimum angle between the fault plane and the raypath plane is 5°, a fully gas hydrate-saturated fault of approximately 70-cm thickness is needed to generate a 0.06-0.08 km/s seismic velocity increase (Figure 10d). This creates inherent nonuniqueness in estimating fault thickness from azimuthal seismic velocity anomalies (Figure 10d). However, we do not see any detectable changes (0.02-0.03 km/s) in seismic velocities from a fault intersecting at a high angle (~30-50°) with raypath planes (for example, the fault shown in green color over azimuthal seismic velocity model for layer LL5 in Figure 6b). We can therefore argue that fault thickness is considerably smaller than 300-400 cm and higher than 20 cm assuming a fully gas hydrate-saturated fault (Figure 10d).

Effect of Faults on Fluid Distribution at Different Scales
In Vestnesa Ridge, the effect of faults on the distribution of gas hydrates has been observed at different scales in different studies. Plaza-Faverola et al. (2015) documented the possibility of a link between tectonic stress and seepage at two 3-D seismic sites separated by approximately 35 km in Vestnesa Ridge. The region with active seepage in Vestnesa Ridge is relatively more faulted compared to other areas where there are pockmarks, but there is no gas seepage. Singhroha et al. (2016) and Singhroha et al. (2019) documented the potential impact of faults on the distribution of gas hydrates in an area covered by 3-D seismic. This paper demonstrates that even at small scales, faults play a big role in spatial distribution of gas hydrates.

Conclusion
This paper documents the application of shooting along circular tracks to study azimuthal seismic velocity variations. We obtain a good and a reliable (with~0.035 km/s uncertainty) azimuthal P-wave velocity model using travel-times from three and four concentric circles in survey1 and survey2, respectively. We find elevated (~0.06-0.08 km/s) seismic velocity anomalies along azimuths in GHSZ where the raypath plane lies along the fault plane or they intersect at a very small angle (<10-20°). High seismic velocity patches occurring along or in the vicinity of inferred fault/fracture systems indicates preferential distribution of gas hydrates along fault/fracture systems in the region. Similarly, we also find reduced (~0.06-0.08 km/s) P-wave velocity anomalies potentially due to the presence of free gas along faults below the BSR.
There are also changes in seismic velocities across faults. Azimuthal variation of seismic velocities document changes in pore-fill as we move across faults and discontinuities. The occurrence of elevated and reduced (>0.05 km/s) P-wave velocities in groups of azimuths bounded by faults indicates compartmentalization of gas hydrates and free gas by fault planes, respectively. Considering the geological setting in Vestnesa Ridge, these results suggest the presence of advection-dominated gas hydrate deposits in faults and the effect of structural features in diffusion-dominated gas hydrate deposits in Vestnesa Ridge. Results from the gas hydrate saturation modeling show that gas hydrate saturated faults of thickness greater than 20 cm and considerably smaller than 300 cm in the layer above the BSR (L5 and LL5) can create observed (0.06-0.08 km/s) azimuthal seismic velocity anomalies. This novel application of azimuthal seismic velocities to study the distribution of gas hydrates shows that the circular shooting can be applied to study variations in seismic velocities at a given site and useful geological information can be derived from it.