Imaging an Underwater Basin and its Resonance Modes using Optical Fiber Distributed Acoustic Sensing

Ambient noise observations using an optical fiber in an underwater basin reveal high-resolution seismic response and resonance effects Shear-wave velocity is imaged via Scholte-wave dispersion and seismic response is revealed via power spectral density and auto-correlation The inferred velocity model is supported by 2D simulations of wave propagation


Introduction
In recent decades, ambient noise tomography has become a common practice in seismology to accurately image subsurface velocity structures (e.g., Shapiro & Campillo, 2004).This objective is typically achieved by extracting Green's functions between pairs of receivers using cross-correlations (CC) of continuous ambient noise recordings (e.g., Foti et al., 2018).These Green's functions, calculated under the assumption of diffuse random wavefields, are then used to resolve the subsurface velocity structure by extracting and analyzing the dispersive properties of surface waves (Park et al., 1999).This method has been successfully applied at global (Nishida et al., 2009) and regional scales (e.g., Guerin et al., 2020;Lin et al. 2008;Mordret et al. 2014).
The scale of ambient noise derived velocity models is currently limited by the geometry of available seismic networks.Seismic sensors are typically deployed at interstation distances of tens to hundreds of kilometers, limiting the spatial resolution of derived velocity models.In addition, underwater seismic recordings are scarce, because of the high cost of deploying and maintaining arrays of ocean-bottom seismometers.Owing to the vast underwater observational gap, only a few studies imaged underwater velocity structures (Guerin et al., 2020;Hable et al., 2019;Mordret et al., 2014;Wolf et al., 2021;Yao et al., 2011).An ideal solution to these two hindrances is the use of distributed acoustic sensing (DAS) on seafloor fiber-optic cables.
The novel technique of DAS enables the transformation of standard optical fibers into dense seismic arrays, continuously recording deformation time-series every few meters across tens of kilometers of fibers.Such measurements have already been utilized for ambient noise tomography both on-land (e.g., Ajo-Franklin et al., 2019;Dou et al., 2017;Yuan et al., 2020) and underwater (Cheng et al., 2021;Spica et al., 2020;Tonegawa et al., 2021;Williams et al., 2021).This high-resolution technique is especially advantageous for ocean-bottom applications since the interrogator is only connected to the on-land end of the fiber, circumventing costly underwater operations.
Unlike previous DAS ocean-bottom ambient noise studies that imaged the shearwave velocity profile of tens of kilometers long sedimentary structures (Spica et al., 2020;Cheng et al., 2021;Tonegawa et al., 2021;Williams et al., 2021), here we apply ambient noise tomography at the scale of few hundreds-of-meters to resolve the complex velocity structure beneath a short 2.5 km-long section of the optical fiber buried in a shallow sedimentary basin offshore Methoni (Greece).This basin displays several rapid lateral variations of the resonance frequency and seismic velocities over a very short distance.Thus, resolving this complex velocity structure is a challenging task that requires the fine spatial sampling of the wavefield achieved by DAS.
To reliably image geological structures, even in the presence of strong lateral variations, we augment the velocity model with power spectral densities (PSD) and autocorrelations (AC) that are useful in understanding local resonance and wave propagation.The PSD is a valuable tool for structural imaging, earthquake hazard mitigation, and engineering purposes as it reveals the basin's response with extremely high resolution and provides a detailed view of frequency-dependent site amplification.Such high-resolution PSDs have been previously generated with waveform simulations in enclosed basins (e.g., Bielak et al., 1999;Sanchez-Sesma et al., 1993;Semblat et al., 2005) but to the best of our knowledge, this study presents the first real-world observations.Similar to the use of CC to obtain the Green's function between different measurement locations, AC represents the 1D Green's function with co-located source and receiver (Claerbout, 1968).Such AC signals have been previously used to map strong impedance contrasts (e.g., Romero & Schimmel, 2018;Taylor et al., 2016), and highlight underwater sedimentary structures and identify fault zones using underwater DAS (Spica et al., 2020;Cheng et al., 2021).
Finally, to examine the inferred shear-wave velocity model, we conducted numerical wave propagation simulations and found good agreement between observations and simulations.The insights obtained from the simulations highlight the strength of underwater DAS in probing wave propagation processes in sedimentary basins, such as scattering, multi-dimensional resonance, and basin edge effects, which are important for the quantification of site effects for seismic hazard assessment.
2 Geological Setting and Fiber Optic Layout and Measurements Methoni basin is located in the very southwest part of Greece and corresponds, for the most part, to the eponymous bay, connected to the Ionian sea.This shallow sedimentary basin is bounded by Methoni village and the mainland in the north and east, Sapientza island in the south, and a natural partially underwater barrier to the west (Figure 1b).The geological settings are described in detail in Text S1.
To image the structure of the Methoni basin, we used underwater DAS ambient noise measurements recorded between April 19th and 25th 2019 by the first 2.5 km of a standard telecommunication optical fiber cable, buried at ~1 m beneath the sediments.As the cable was buried in shallow waters, its route is relatively well-constrained.The measurement setup is detailed in text S2.
When entering the sea, the cable is deployed over Alluvium deposits (Fytrolakis, 1980) (Figure 1b), and it is buried under soft underwater sediments until it reaches the Limestone formations at the western edge of the basin.The extension of the Flysch formation, which covers the eastern edge of the basin and Nisakouli island (Figure 1b) (Fytrolakis, 1980), onto the seabed is unknown, and it may exist beneath limited cable sections.It is also expected that the Limestone formations extend beneath the Alluvium and Flysch, and result in a sharp velocity contrast.

Extracting Scholte-Wave Dispersion curves
DAS recordings in Methoni bay are dominated by surface gravity waves and ambient noise Scholte waves.The latter is used here to image the structure beneath Methoni bay.These Scholte-waves are observed here at frequencies between 1 and 5 Hz (e.g.Sladen et al., 2019), and are locally amplified and modulated by the sedimentary structure beneath the fiber (Lior et al., 2021a).
To estimate the velocity structure, we retrieve empirical Green's functions by calculating CC on short cable sections and pick Scholte-wave phase-velocity dispersion curves in the frequency-wavenumber (FK) domain.To resolve small-scale lateral structure variations while maintaining reasonable wavenumber resolution, short cable segments of 31 and 15 equally spaced channels (595.2 m and 288 m, respectively) were used to compute CC both in oceanward and landward directions.Green's function retrieval and dispersion curve picking are detailed in text S3.For many locations along the fiber, each of the four CC gathers and FK transforms reveals slightly different dispersion characteristics.Such differences are shown in Figure 2. The shoreline is located at channel 19 (~365 m), so part of the longer (31 channels) segment is located on-land.Landward CC gathers show reflections from the edge of the basin (Figure 2e and 2g from 3 s), suggesting a sharp impedance contrast.Similar behavior is observed at the other edge of the basin (not shown).Furthermore, the energy of the longer oceanward segment (Figure 2a and 2b) is lower compared to the other CC gathers since the virtual source is located on-land.In addition, for landward FK (Figure 2f and 2h), the first higher mode is significantly more visible compared with oceanward FK (Figure 2b and 2d).Thus, considering all four gathers when picking dispersion curves provides a more complete image of Scholte-wave propagation, as indicated by the blue and black circles in Figure 2 (panels b, d, f and h), which correspond to the picked fundamental and first higher modes, respectively.

Retrieving a 2D Velocity Model from Scholte-wave dispersion curves
The 2D shear-wave velocity profile beneath the fiber is resolved in two stages.First, we invert for a one-dimensional (1D) profile beneath each short segment of the fiber.Then, we select one location where a reliable profile was obtained and invert for its neighboring profiles, allowing for small changes (10%) with respect to adjacent, previously estimated profiles.The second step is performed recursively to obtain a smoothed 2D velocity model.The inversion algorithm is based on the Neighbourhood Algorithm (NA) developed by Sambridge (1999), and the implementation introduced by Mordret et al. (2014), as further described in Text S4.
Each 1D shear-wave velocity model is composed of a water layer, three solid layers and a bottom half-space; Overall we invert for seven parameters: three layer thicknesses and four velocities.We found that this configuration provides sufficient spatial resolution in depth while preventing model overfitting.As the starting point for the 2D model, we chose channel 32 (615 m) along the fiber, previously displayed in Figures 2, where the fundamental and first higher modes are resolved at a broad frequency range, and nearly identical velocity profiles were obtained for ten NA inversions initialized with different starting models (red curves in panels b, d, f and h of Figure 2 and Figure S1).Since the velocity of the bottom half-space is poorly constrained (1617 ± 206 / ), we turn to the on-land seismic stations (Figure 1b) located on the limestone formations, assuming that they are located on the same geological formation found at the bottom of the basin.This assumption is reasonable given the geological setting described in Section 2. CC calculations on ambient noise recorded by the two on-land stations (Figure 1b) are described in Text S4.This analysis suggests a shear-wave velocity of 1811 ± 73 / , consistent with that expected for fractured limestone.
To obtain a 2D velocity profile, the inversion was recursively applied along both fiber directions (oceanward and landward), starting at channel 32, with fixed bottom halfspace velocity of 1.8 km/s, as estimated from on-land CC.The smoothed 2D velocity model is plotted in Figure 3a.The model displays several prominent features: 1) a steep velocity contrast between the bedrock (Vs=1800 m/s) and the overlying layer (Vs=468 ± 90 / ), 2) two deeper sediment regions around channels 50 and 80, and 3) sedimentary layers are generally thinner towards the basin edges.

Power Spectral Density and Auto-Correlation
To further investigate the basin geometry, strain-rate PSDs and AC were calculated as detailed in Text S5.The PSD

Earthquake Observations
To support the observed AC (Figure 3c), we investigate wave propagation characteristics of earthquakes recorded by the fiber.Since earthquake recordings are a convolution of the source time-function and the Green's function between source and receiver positions, near-vertical seismic arrivals are closely related to the 1D Green's functions displayed in the AC image (Figure 3c).We examined four different local earthquakes that occurred at different locations and source-fiber azimuths (Lior et al., 2021a) for which near-vertical direct body wave arrivals are expected.Apart from SNR variations, wave propagation is extremely similar for the different earthquakes, suggesting that local effects dominate these recordings.The first five seconds following the direct Swave arrival of a magnitude 2 at a hypocentral distance of 37 km are plotted in Figure 3d.These signals significantly resemble the AC in Figure 3c.

Numerical simulations
To support the obtained velocity model (Figure 3a), we conducted 2D numerical wave propagation simulations in a simple, smooth basin model, as previously done by Lior et al. (2021b).These simulations cannot capture the full propagation characteristics expected from 3D wave propagation in such a complex real medium, recorded by a curved sensor array (Figure 1b); They were conducted to qualitatively recreate the general features resolved by the velocity model in Figure 3 and understand their relation to wave propagation.In particular, the patterns of seismic response of a basin, observed here via PSD and AC analysis, results from both lateral and vertical resonances (e.g., Bielak et al 1999), which are partially captured by a 2D simulation.
Two velocity models are tested (Figures S4a and 4a), which are smoothed versions of that shown in Figure 3a.The two models differ in the manner in which the depth of the bottom layer varies surrounding channel 80. Since the obtained 2D model (Figure 3a) is affected by the inversion approach, which allows for small model variations between adjacent 1D profiles (Section 3.1.2),the estimated layer interface may in fact be steeper.Thus, in Figure 4a we model a steeper increase in depth with respect to that in Figure S4a.In these modified models, each layer has uniform seismic velocities, set as the average of the layer velocities in the original model (Figure 3a), and smoothed layer interfaces.The velocity model in Figure 3a does not show obvious basin edges at the two ends of the model, most likely due to the rapid lateral variations and smoothing constraints of the 2D inversion.However, observed PSD and AC show clear basin edge patterns identified as increasing frequencies (Figure 3b) and scattered waves (Figure 3c), respectively.To qualitatively reproduce the observed basin edge patterns, the models in Figures S4a and 4a were manually closed on both lateral ends to generate strong impedance contrasts at the basin's edges that enable lateral resonance in the basin.The simulation set-up is further described in Text S6.S4b and 4b) and AC (Figures S4c and 4c) are obtained as described in Text S5.To highlight the same frequencies as in the observed PSD (Figure 3b), the simulated PSDs were divided by the spectrum of the Ricker wavelet used as source time function in the simulations (see Text S6) and multiplied by the average (along the spatial axis) observed (Figure 3b) PSDs.For comparison with the earthquake recording in Figure 3d, the first five seconds of the simulation are shown in Figures S4d and 4d

Discussion
Each individually-inverted 1D profile represents the velocity beneath its associated CC gather, assuming negligible lateral variations.Thus, each 1D profile should be treated with caution, as it may be biased by lateral variations within the CC fiber segment (at least 15 channels, i.e., 288 m).For fiber segments with strong lateral variations, the information carried by the PSD and AC is crucial, as discussed further.
Several features in the estimated velocity model (Figure 3a) have clear signatures in the observed PSD (Figure 3b), AC (Figure 3c) and earthquake recording (Figure 3d).Sedimentary layers generally get thinner towards the edges of the basin (panel a), causing reflections seen in both the AC (panels c) and earthquake record (panel d).These reflected waves propagate at velocities of 400 to 600 m/s, as indicated by the reference 500 m/s red curves.Basin edge effects, a result of interference between vertical and lateral resonances, also appear in the PSD (panel b) as a frequency increase towards the basin's edges, as predicted by wave-propagation simulation studies (e.g., Papageorgiou & Kim, 1991;Sanchez-Sesma et al., 1993).The deeper sediment region around channel 80 (panel a) is associated with a striking triangular shape in the AC and earthquake recording, as well as a local PSD frequency increase, suggesting scattered waves.Another clear feature in the AC and earthquake signal are late (7 to 10 seconds) landward propagating waves, seen from channel 40 to 80 (panel c).These cannot be correlated with any feature in the velocity model (panel a), and may be a result of a deeper velocity structure.
At the basin's edges, significant lateral variations and a sharp impedance contrast result in horizontally reflected waves (e.g., CC gathers in Figure 2) as well as abrupt disruption of seismic signals at channels 19 (shoreline) and 145 (Figure 3b-3d).The latter may also be a result of unfavorable cable-ground coupling.This impedance contrast is not properly resolved in the 2D profile (Figure 3a).To address this issue and improve the physical reliability of the simulations, the models were manually closed ("pinned-up") on both sides (Figures S4a and 4a), and reliably reproduce the observed reflections (Figure S4c-d and 4c-d).It should be noted that simulated waves (Figure S4d and 4d) reflected from the left (on-land) edge of the basin exhibit dispersive characteristics, not observed in the earthquake recording (Figure 3d).
The waveform simulations in Figures S4 and 4 are found to be in good agreement with observations (Figure 3), supporting the estimated 2D velocity model.Specifically, the simulations that exhibit a steeper depth change for the bottom layer around channel 80 (Figure 4) better reproduced the observations.This agreement is observed for waves reflected from the basin edges, where observations and simulations show similar velocities (Figure 3c and 4c, and 3d and 4d).The increasing frequencies at the basin edges are also well reproduced in the PSD.The AC and PSD features seen around channel 80 (Figure 3b-3c) are qualitatively reproduced by the simulations (Figure 4b-4c), a result of complex wave propagation associated with the sharp increase in the depth of the bottom layer.As expected, waves seen at 7 to 10 seconds in Figure 3c (channels 40 to 80) are not reproduced in Figure 4c, further suggesting that these are related to a deeper velocity structure.
The PSD and AC features observed around channel 80 (Figure 3b-3c), are caused by lateral wave propagation effects that may be associated with several different structural phenomena.This structure may be an underwater canyon, filled with low-velocity sediments, that generates reflected waves from its edges.Another option is hinted by the proximity of channel 80 to Nisakouli island (Figure 1), where a Flysch formation is outcropping on the island (Fytrolakis, 1980).The Flysch may extend underwater beneath the fiber and appear as a low-velocity outcrop overlaying the limestone formation (bottom layer).The relation between the different formations may play a role in generating 3D wave propagation complexities associated with the observations.A third possibility is that this region corresponds to a localized fault zone, characterized by low velocity fractured rocks.Such a structure would display trapped high-frequency waves and generate an indicative triangular wave propagation pattern (e.g., Lindsey et al., 2019;;Spica et al., 2020).Also, an increase in the PSD frequency content, as seen here (Figure 3b), can originate from wave propagation effects at the boundary between two sub-basins (e.g., Bielak et al 1999), but may equivalently be generated if a curved fiber approaches any edge of the basin.Such boundaries may appear in the velocity model (Figures 3a, S4a and 4a) at channels ~60 and ~90 for the left and right sub-basins, respectively, yet these are too subtle and cannot reproduce the AC and PSD observations (Figure 3b-3c) as demonstrated by the simulations in Figure S4.With a significant lack of ground truth, we cannot determine which of these hypotheses is the valid one.Further simulations and broader scale imagery should be conducted to accurately resolve this structure and understand local wave propagation.

Conclusions
We demonstrated the ability to image the shear-wave velocity model beneath a short underwater optical fiber using continuous ambient noise recordings.AC and PSD provide a detailed view of wave propagation effects, and basin resonance and amplification patterns, respectively.To the best of our knowledge, this is the first time that such high-resolution AC and PSD were obtained using real-world observations.A smoothed 2D velocity model was obtained by analyzing the dispersion of Scholte-waves using short fiber segments.This model was modified following the insights gained from the PSD and AC analysis, and subsequent waveform simulations revealed a good agreement with observations.However, additional knowledge and constraints on local underwater geology are required to reliably interpret the obtained model.
The imagery approach used here is the first step towards a method that combines surface wave dispersion as well as AC and PSD analysis using DAS fiber optic data.Such a method will account for both lateral and vertical resonance to reliably and objectively image complex structures in 2D and allow for greater sensitivity to strong lateral variations.For a more complete 3D model, additional broadband seismometers and/or augmenting optical fibers should be used.Waveform simulations proved to be an excellent tool to validate the inferred model, and should be used and expanded to 3D as needed.

Figure 1 .
Figure 1.Map of the study area.(a) South Greece, study area is marked by red lines.(b) Aerial map of Methoni basin, the fiber location is color coded by water depth and channel numbers and distances along the fiber are indicated.On-land seismometers are marked by black rectangles.

Figure 2 .
Figure 2. CC gathers and corresponding FK transforms around channel 32 calculated: oceanward using (a-b) 31 and (c-d) 15 channel fiber segments, and landward using (e-f) 31 and (gh) 15 channel fiber segments.FK transforms in panels (b, d, f and h) are calculated for the CC gathers in panels (a, c, e and g) following zero-padded.The picked fundamental and first higher modes are indicated by blue and black circles, respectively (panels b, d, f and h).Theoretical dispersion curves associated with the best fitting shear-wave profiles are plotted in red curves for each of the ten different inversions (panels b, d, f and h).These represent the velocity model beneath channel 32 along the fiber.The colorbar corresponds to panels (a, c, e and g).
Figure 3. (a) 2D shear wave velocity profile along the fiber, (b) Amplitude corrected PSD, (c) AC following FK filter to remove velocities lower than 800 m/s (d) 5 seconds following the direct Swave of a magnitude 2 earthquake recorded at a distance of 37 km, bandpass filtered between 1 and 3 Hz.In (a), the velocity of the bottom half space and water layer are fixed at 1.8 km/s and 0, respectively.Red lines in (c) and (d) indicate velocities of 500 m/s.The shoreline is indicated by black vertical lines in (a)-(c). .

Figure 4 .
Figure 4. Simulation results for a steep increase of bedrock depth around channel 80.(a) Smoothed 2D shear-wave velocity profile.Layer velocities are indicated in the legend.(b) PSD and (c) AC calculations.PSD are convoluted with the average observed PSD (see text).An FK filter is applied to the AC in (c) to remove velocities lower than 800 m/s.(d) The first 5 seconds of the simulated waveforms.Red lines in (c) and (d) indicate velocities of 500 m/s.