Loïc Viens
email: lviens@umich.edu

Luis Fabian Bonilla

Zack J Spica

Kiwamu Nishida

Tomoaki Yamada

Masanao Shinohara

Nonlinear earthquake response of marine sediments with distributed acoustic sensing

AutoCorrelation Functions (ACFs) of earthquakes recorded by a DAS experiment 10 exhibit phase delays with increasing ground motions 11 • ACF phase delays are converted to relative velocity drops in the medium that char-12 acterize soil nonlinearity 13 • DAS is used to infer the nonlinear behavior of soils with an unprecedented spa-14

Introduction

Local geological conditions can significantly impact the propagation of incoming seismic waves from earthquakes. In particular, shallow, soft, and unconsolidated sediment layers are well known to amplify earthquake ground motions [START_REF] Sanchez-Sesma | Site effects on strong ground motion[END_REF], which can lead to catastrophic events such as during the 1985 moment magnitude (M w ) 8.0 Michoacán earthquake in Mexico [START_REF] Anderson | Strong ground motion from the Michoacan, Mexico, earthquake[END_REF][START_REF] Campillo | Destructive strong ground motion in Mexico city: Source, path, and site effects during great 1985 Michoacán earthquake[END_REF]. When subjected to weak dynamic strains (i.e., less than 10 -4 and 10 -8 for field observations and laboratory experiments, respectively; [START_REF] Ishihara | Soil behaviour in earthquake geotechnics[END_REF][START_REF] Tencate | Nonlinear and nonequilibrium dynamics in geomaterials[END_REF], shallow soils linearly amplify seismic waves. During large dynamic strains, however, soft sediments can behave nonlinearly (e.g., [START_REF] Field | Nonlinear ground-motion amplification by sediments during the 1994 northridge earthquake[END_REF][START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF]. Soil nonlinearity is generally characterized by a relative reduction of the high-frequency groundmotion amplification, which is related to an increase of damping in the medium, and a shift of the resonance frequency to lower frequencies due to a reduction of the shear modulus [START_REF] Beresnev | Nonlinear soil response-a reality?[END_REF][START_REF] Brunet | Transitional nonlinear elastic behaviour in dense granular media[END_REF][START_REF] Bonilla | Nonlinear site response evidence of K-NET and KiK-net records from the 2011 off the Pacific coast of Tohoku Earthquake[END_REF][START_REF] Lyakhovsky | Non-linear damage rheology and wave resonance in rocks[END_REF][START_REF] Zaitsev | Luxemburg-gorky effect in a granular medium: Probing perturbations of the material state via cross-modulation of elastic waves[END_REF]. In some cases, large dynamic strains can trigger a complete failure of cohesionless and saturated shallow sediments through soil liquefaction [START_REF] Kramer | Geotechnical Earthquake Engineering[END_REF], which can have disastrous consequences for human infrastructures as observed during the 1964 Niigata (Japan, [START_REF] Ohsaki | Niigata earthquakes, 1964 building damage and soil condition[END_REF][START_REF] Beyreuther | ObsPy: A Python Toolbox for Seismology[END_REF]-2011Christchurch (New Zealand, Quigley et al., 2013) earthquakes. Therefore, characterizing the nonlinear response of shallow sediments to earthquake ground motions is critical for mitigating seismic risk.

Several empirical methods have been developed to assess the response of soils to ground motions. A classical approach relies on computing the spectral ratio of earthquakes recorded at a soft-soil site and a nearby reference rock site [START_REF] Borcherdt | Effects of local geology on ground motion near san francisco bay[END_REF]Field & Ja-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Geophysical Research Letters cob, 1995; [START_REF] Bonilla | Site amplification in the san fernando valley, california: Variability of site-effect estimation using the s-wave, coda, and h/v methods[END_REF]. However, a reference site may not always be available in the vicinity of the site of interest. Another approach consists in using pairs of surfaceborehole stations to detect potential soil nonlinear elastic behavior between the two sensors [START_REF] Bonilla | Nonlinear site response evidence of K-NET and KiK-net records from the 2011 off the Pacific coast of Tohoku Earthquake[END_REF][START_REF] Minato | Monitoring seismic velocity change caused by the 2011 tohoku-oki earthquake using ambient noise records[END_REF][START_REF] Nakata | Near-surface weakening in japan after the 2011 tohoku-oki earthquake[END_REF][START_REF] Régnier | Assessing nonlinear behavior of soils in seismic site response: Statistical analysis on kiknet strongmotion data assessing nonlinear behavior of soils in seismic site response[END_REF]Sawazaki et al., 2006;[START_REF] Takagi | Coseismic velocity change in and around the focal region of the 2008 Iwate-Miyagi Nairiku earthquake[END_REF][START_REF] Wen | Investigation of non-linear site amplification at two downhole strong ground motion arrays in taiwan[END_REF]. While this technique isolates the shallow subsurface response from the earthquake source and path effects, pairs of surface-borehole instruments are expensive to install and their low spatial coverage does not capture small-scale lateral variations.

AutoCorrelation Functions (ACFs) calculated from data recorded by surface seismometers yield the reflectivity response of the underlying elastic structure [START_REF] Claerbout | Synthesis of a layered medium from its acoustic transmission response[END_REF][START_REF] Wapenaar | Synthesis of an inhomogeneous medium from its acoustic transmission response[END_REF]. This technique has been primarily applied to image interfaces with strong seismic impedance contrasts using earthquake [START_REF] Delph | Constraining crustal properties using receiver functions and the autocorrelation of earthquake-generated body waves[END_REF][START_REF] Pham | On the feasibility and use of teleseismic p wave coda autocorrelation for mapping shallow seismic discontinuities[END_REF][START_REF] Tork Qashqai | Crustal imaging with Bayesian inversion of teleseismic P wave coda autocorrelation[END_REF][START_REF] Viens | Imaging the Kanto Basin seismic basement with earthquake and noise autocorrelation functions[END_REF] and ambient seismic field (ASF; [START_REF] Gorbatov | Crustal properties from seismic station autocorrelograms[END_REF][START_REF] Ito | Autocorrelation analysis of ambient noise in northeastern Japan subduction zone[END_REF][START_REF] Kennett | Lithosphere-asthenosphere p-wave reflectivity across australia[END_REF][START_REF] Saygin | Retrieval of the P wave reflectivity response from autocorrelation of seismic noise: Jakarta Basin, Indonesia[END_REF][START_REF] Spica | Marine sediment characterized by ocean-bottom fiber-optic seismology[END_REF][START_REF] Viens | Imaging the Kanto Basin seismic basement with earthquake and noise autocorrelation functions[END_REF] datasets recorded onshore.

Repeated ACF computations through time from continuous ASF time series have also been used to monitor temporal seismic velocity changes in the subsurface in different environments, such as volcanic [START_REF] De Plaen | Single-station monitoring of volcanoes using seismic ambient noise[END_REF][START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at merapi volcano, indonesia[END_REF][START_REF] Yates | Volcanic, coseismic, and seasonal changes detected at white island (whakaari) volcano, new zealand, using seismic ambient noise[END_REF] and earthquake source [START_REF] Hobiger | Single-station crosscorrelation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 iwate-miyagi nairiku earthquake[END_REF][START_REF] Ohmi | Temporal variations of crustal structure in the source region of the 2007 noto hanto earthquake, central japan, with passive image interferometry[END_REF][START_REF] Wegler | Sudden drop of seismic velocity after the 2004 mw 6.6 mid-niigata earthquake, japan, observed with passive image interferometry[END_REF] regions. However, the partitioning of body and surface waves in ACFs computed from the ASF is generally unknown and hinders the interpretation of the measured velocity changes [START_REF] Nakahara | Auto correlation analysis of coda waves from local earth[END_REF]. To ease the interpretation, ACFs have also been computed from earthquake P-, S-, or coda-wave windows [START_REF] Bonilla | Monitoring Coseismic Temporal Changes of Shallow Material during Strong Ground Motion with Interferometry and Autocorrelation[END_REF]Bonilla & Ben-Zion, 2020;[START_REF] Nakahara | Auto correlation analysis of coda waves from local earth[END_REF][START_REF] Qin | Imaging and monitoring temporal changes of shallow seismic velocities at the garner valley near anza, california, following the m7.2 2010 el mayor-cucapah earthquake[END_REF]. [START_REF] Bonilla | Detailed space-time variations of the seismic response of the shallow crust to small earthquakes from analysis of dense array data[END_REF] showed that the first negative peak of ACFs calculated during earthquake ground motions corresponds to the seismic-wave two-way travel time between the sensor and the first major interface below the station, and captures the soil nonlinear response. Moreover, [START_REF] Bonilla | Monitoring Coseismic Temporal Changes of Shallow Material during Strong Ground Motion with Interferometry and Autocorrelation[END_REF] and [START_REF] Qin | Imaging and monitoring temporal changes of shallow seismic velocities at the garner valley near anza, california, following the m7.2 2010 el mayor-cucapah earthquake[END_REF] showed that the response of the shallow subsurface obtained from ACFs at surface stations yields a similar estimation of the soil nonlinear behavior as that from a surface-borehole station configuration. In other words, ACFs can isolate the site response term from the earthquake source and path effects, making single-component stations a powerful tool for analyzing shallow sediment nonlinear behavior.

Mapping local site effects with data-driven techniques remains challenging due to the large density of seismometers needed to capture complex spatial variations of the seismic wavefield. In some cases, a large station coverage can be nearly impossible to attain due to environmental or physical constraints, especially in urban and underwater areas.

Nevertheless, recent technological advances in Distributed Acoustic Sensing (DAS) offer an unprecedented opportunity to measure the Earth's vibrations over tens of kilometers with a dense spatial resolution (∼1-50 m) by turning ground-coupled fiber-optic cables into arrays of sensors [START_REF] Hartog | An introduction to distributed optical fibre sensors[END_REF]. DAS uses an optoelectrical interrogator to probe fibers with a laser sending thousands of short pulses of light every second. As each pulse of light travels down the fiber, some of the light is reflected back to the interrogator in a process known as Rayleigh backscattering. External forcing, such as seismic waves, generate phase shifts of the back-scattered Rayleigh light, which are measured by the interrogator. The measured phase shifts are finally linearly converted to longitudinal strain (or strain rate) along the cable over a sliding spatial distance (i.e., the gauge length). Both fit-to-purpose and existing fiber-optic telecommunication cables have been used to record high-fidelity earthquake wavefields [START_REF] Lellouch | Seismic velocity estimation using passive downhole distributed acoustic sensing records: Examples from the san andreas fault observatory at depth[END_REF][START_REF] Lior | On the Detection Capabilities of Underwater Distributed Acoustic Sensing[END_REF][START_REF] Spica | Urban seismic site characterization by fiber-optic seismology[END_REF][START_REF] Spica | Subsurface imaging with ocean-bottom distributed acoustic sensing and water phases reverberations[END_REF][START_REF] Wang | Ground motion response to an ML 4.3 earthquake using co-located distributed acoustic sensing and seismometer arrays[END_REF][START_REF] Zeng | Properties of Noise CrossCorrelation Functions Obtained from[END_REF]. One great advantage of telecommunication fibers is that they have been widely deployed, from the oceans' bottom to nearly every street in large developed cities, to sustain our modern telecommunication network. Therefore, DAS could complement expensive urban and off-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Geophysical Research Letters shore seismic array deployments by probing existing telecommunication cables to capture the full extent of earthquake wavefields.

Nonlinear behavior of marine sediments during earthquake ground motions has been observed with sparse ocean bottom seismometers [START_REF] Dhakal | Assessment of nonlinear site response at ocean bottom seismograph sites based on s-wave horizontal-to-vertical spectral ratios: a study at the sagami bay area k-net sites in japan[END_REF][START_REF] Kubo | GroundMotion Characteristics and Nonlinear Soil Response Observed by DONET1 Seafloor Observation Network during the 2016 Southeast OffMie, Japan, Earthquake[END_REF].

This study focuses on capturing the fine spatial resolution of the response of shallow marine sediments to 105 earthquakes recorded along a telecommunication cable offshore the Sanriku coast in Japan by a DAS experiment (Figures 1a-b). After filtering the data into different frequency bands, we calculate ACFs from the earthquake ground motions to infer the soil response at different depths. We detect changes in the ACF time series that are converted to relative velocity changes to characterize the soil linear and nonlinear regimes below each DAS channel. We show that the relative velocity changes exhibit spatial variations along the cable and are frequency dependent. We finally discuss our results and the potential of DAS for extracting soil parameters with unprecedented spatial resolution.

2 Data and Methods

DAS data

The Earthquake Research Institute, The University of Tokyo, operates an oceanbottom observatory composed of three 3-component accelerometers and two tsunami meters offshore the Sanriku Coast (Figure 1b; [START_REF] Kanazawa | Ocean-bottom observatory for earthquakes and tsunami off sanriku, north-east japan using submarine cable[END_REF]Shinohara et al., 2021Shinohara et al., , 2022)). The data recorded by the instruments are streamed in real-time to the landing station located in the city of Kamaishi through a submarine telecommunication cable. The cable contains six dark (unused) dispersion-shifted single-mode optical fibers with a wavelength of 1,550 nm, which are suitable for DAS measurements. The first 47.7 km of the cable are relatively straight and buried under 0.6-0.7 m of sediments, which should guarantee a good coupling with the surrounding medium. Moreover, the relatively weak levels of shaking analyzed in this study (i.e., peak strain < 10 -7 or peak ground acceleration < 1 cm/s/s) should not generate fiber movements that could impact the recorded seismic waves [START_REF] Nakamura | Rotation motions of cabled oceanbottom seismic stations during the 2011 Tohoku earthquake and their effects on magnitude estimation for early warnings[END_REF]).

An AP Sensing N5200A DAS interrogator unit [START_REF] Cedilnik | Ultralong reach fiber distributed acoustic sensing for power cable monitoring[END_REF] probed one of the dark fibers between November 18 and December 2, 2019, and recorded continuous data over the first 70 km of the cable with a sampling rate of 500 Hz. The gauge length and spatial sampling are set to 40 m and 5.1 m, respectively. The DAS system recorded hundreds of earthquakes during the two weeks of measurement. We first convert the raw DAS data to strain [START_REF] Shinohara | Performance of seismic observation by distributed acoustic sensing technology using a seafloor cable off sanriku, japan[END_REF] and focus on the ground motions from 105 earthquakes that are clearly recorded by all DAS channels (Figure 1a-b). The velocity magnitude (M V , determined from the amplitude data recorded by short-period velocity seismometers; [START_REF] Kanbayashi | A method for determining magnitude of shallow earthquakes occurring in and near japan[END_REF] of the earthquakes ranges between 1.0 and 5.3, and we show the strain waveforms of a M V 2.5 earthquake in Figure 1c. This event occurred on November 28, 2019 at 14:17:32UTC at a depth of 30 km. Clear P-and S-wave arrivals can be observed at most channels as well as locally generated surface waves that significantly extend the ground motion duration.

Autocorrelation functions and relative velocity changes

For each earthquake and each DAS channel, we compute the time derivative of the strain data to retrieve strain-rate waveforms, which are proportional to acceleration time series. We bandpass filter the strain-rate data between 2 and 30 Hz (all filters are twopass four-pole Butterworth bandpass filters) and select a fixed 15-s window starting 5 s before the earthquake absolute maximum amplitude. We then bandpass filter the raw strain-rate waveforms into 18 frequency bands (e.g., [3][4][5][6][4][5][6][7][8]...,, apply 1-bit normalization to the data [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[END_REF], and compute ACFs over the fixed 15-

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Geophysical Research Letters s window using a phase correlation method in the frequency domain [START_REF] Schimmel | Noise reduction and detection of weak, coherent signals through phase-weighted stacks[END_REF][START_REF] Ventosa | Towards the Processing of Large Data Volumes with Phase CrossCorrelation[END_REF]. We show the bandpass filtered strain-rate waveforms of the M V 2.5 earthquake together with their corresponding ACFs at channel 5000 in Figures 1d and1e, respectively.

As the ACFs are calculated around the S-wave direct arrival, we expect their first negative peak to capture the S-wave two-way travel time between the cable and an interface with a strong seismic impedance contrast (Bonilla & Ben-Zion, 2020). Nevertheless, a theoretical/numerical study should explore the potential effects of ACF near-zerotime-lag source effects on the S-wave two-way travel time, especially in offshore environments. In addition, the different frequency bands allow us to sample different depth of the media, with low-frequency bandpass filtered ACFs displaying later arrivals as they sample deeper media compared to high-frequency ACFs. However, detailed knowledge of the subsurface is required to precisely determine the depth sampled by ACFs.

The 105 earthquake waveforms analyzed in this study generated various levels of dynamic strain along the cable. In Figures 2a-d, we show the ACFs calculated for all the earthquakes after bandpass filtering the strain-rate data between 10-20 Hz and 15-30 Hz at channels 5000 and 7000. We also show the dynamic peak strains computed as the absolute maximum amplitude of the bandpass filtered strain data in Figures 2e-f. For both frequency bands, the ACF first negative peaks exhibit similar lag times for weak dynamic strains (<∼5×10 -10 ), but clear delays can be observed for larger dynamic peak strains.

Soil nonlinear behavior during ground motions delays the ACF first negative peak and can therefore be interpreted as a velocity reduction in the medium (Bonilla & Ben-Zion, 2020). Under the assumption that the changes in the medium are uniformly distributed, we can estimate the relative velocity changes (dv/v) of each ACF with respect to a reference ACF with the stretching method [START_REF] Lobkis | Coda-wave interferometry in finite solids: Recovery of p-to-s conversion rates in an elastodynamic billiard[END_REF][START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at merapi volcano, indonesia[END_REF] as

τ = dt t = - dv v , (1)
where τ , dt/t, and dv/v are the stretching coefficient, the relative time shift, and the relative velocity change, respectively. The dv/v changes are computed in several steps to estimate measurement errors. For each channel and each frequency band, we first compute 100 reference ACFs by randomly stacking between 5 and 20 ACFs computed during the 20 weakest dynamic peak strains. We then select a time window that corresponds to 75% of the inverse of each lower cutoff frequency (e.g., the first 0.25 s of the ACF for the 3-6 Hz frequency band) to focus on the first negative peak of the ACFs. We stretch/compress the selected window of each reference ACF to find the stretching coefficient that maximizes the fit with each earthquake ACF, and therefore measure dv/v changes. The stretching is performed in two steps; we first use ten values uniformly distributed between -50 and 50% of stretching to find an initial guess of the stretching coefficient, and then refine the measurement by interpolating the stretched waveforms 500 times between the neighboring values (similar to [START_REF] Viens | Complex near-surface rheology inferred from the response of greater tokyo to strong ground motions[END_REF]. For each channel and each frequency band, we finally compute the mean and one standard deviation from the mean of the 100 dv/v measurements to estimate errors.

Results

We show the mean of the dv/v measurements for all the frequency bands and earthquakes for two ranges of channels in Figures 3a-b, and the corresponding uncertainties in Figure S1. While the soil nonlinear response can rapidly evolve spatially, we display the combined results at 10 neighboring channels (i.e., over 51 m) between channels 5000-5010 and 7000-7010 for visibility. For dynamic peak strains smaller than ∼5×10 -10 , dv/v measurements are generally equal to zero for all frequency bands at both locations, which

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Geophysical Research Letters indicates that there is no change in the medium. However, clear dv/v drops can be observed in different frequency bands with increasing dynamic peak strains at the two locations. For example, we primarily observe dv/v reductions between central frequencies (i.e., the central frequency of the bandpass filter: 15 Hz for the 10-20 Hz bandpass filter) of 15-24 Hz for channels 5000-5010 and 12-28 Hz for channels 7000-7010. Moreover, we also note that the intensity of the dv/v changes varies, with larger changes at channels 5000-5010 compared to those at channels 7000-7010.

Spatial variations of the relative velocity changes can also be detected along the cable with the high density of DAS channels. In Figures 3c-d, we show the mean dv/v measurements along the cable in two frequency bands (5-10 and 10-20 Hz), which exhibit clear spatial differences. In the 5-10 Hz frequency band, almost no dv/v changes can be observed between channels 500 and 6900, even during the strongest dynamic peak strains. However, we detect clear dv/v drops for dynamic peak strains above 10 -9 between channels 6900 and 8200. In the 10-20 Hz frequency band, almost no changes are found between channels 500-3500, but large dv/v drops are observed after channel 3500 for dynamic peak strains larger than ∼ 10 -9 .

To isolate and investigate the average sediment response during strong ground motions, we also compute dv/v measurements between each reference ACF (i.e., the stack of the ACFs computed during the 20 weakest dynamic peak strains) and a stack of the ACFs computed during the 10 largest dynamic peak strains (More details about uncertainty computation in Text S2). Mean dv/v changes between weak and strong ground motion ACFs exhibit clear spatial and frequency variations (Figure 4a). Between channels 500 and 2000, we do not detect any large dv/v changes. However, we observe spatial variations of the dv/v reductions at central frequencies above 15 Hz between channels 2000 and 9000. We also observe clear dv/v changes at frequencies below 15 Hz between channels 6300 and 9000. Note that the largest dv/v values are not correlated with the highest dynamic peak strain values (Figure 4b). Therefore, the coherent spatial changes across frequency bands highlight the sensitivity of DAS ACFs to local site conditions.

The amplitude of dv/v reductions is expected to increase with increasing dynamic peak strains. In Figures 4c-e, we show the dv/v measurements calculated between the weak and strong ground motion ACFs as a function of the dynamic peak strains in three frequency bands (3)(4)(5)(6)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30). We only show the results at 650 locations between channels 1700-8200 as we average the dv/v and dynamic peak strain values over 10 neighboring channels (e.g., channels 1995-2005 for channel 2000). This step is performed to compare our results with local site condition data from a velocity model of the region, as discussed below. For the three frequency bands, the largest dynamic peak strains are typically observed before channel 3000 and we generally see an increase of the dv/v drops with increasing dynamic peak strains. In the 15-30 Hz frequency band, the dv/v values at the beginning of the cable (before channel 4000) do not behave in the same way as channels located further away from the coast, as shown by the relatively scattered dv/v values in Figure 4d.

In Figures 4f-h, we show the dv/v measurements at the same 650 locations as a function of the average S-wave velocity in the first 30 m of the ground (V S30 ) obtained from the velocity model derived by [START_REF] Viens | Understanding surface-wave modal content for high-resolution imaging of submarine sediments with distributed acoustic sensing[END_REF]. The [START_REF] Viens | Understanding surface-wave modal content for high-resolution imaging of submarine sediments with distributed acoustic sensing[END_REF] model was obtained by inverting Rayleigh wave phase velocity dispersion curves calculated by seismic interferometry using virtual sources located every 10 channels (e.g., 51 m). We compute V S30 from the 650 locations of the velocity model and apply a smoothing of the V S30 values over 5 locations. We observe a decrease of dv/v changes with decreasing V S30 values in the 20-40 Hz frequency band. However, we do not observe any clear correlation between V S30 and the dv/v results in the 3-6 and 15-30 Hz frequency bands.

Accepted Article

This article is protected by copyright. All rights reserved. While larger dynamic peak strains generally correlate with larger dv/v drops, the correlation with V S30 is weaker or even nonexistent. This contrasts with the study by [START_REF] Yang | Sub-kilometer correlation between near-surface structure and ground motion measured with distributed acoustic sensing[END_REF], who observed a good correlation between stronger site effects and slower shallow S-wave velocities. Three hypotheses can explain this behavior. First, the velocity model, which was obtained from ASF cross-correlation functions spanning over 2 km (i.e., 400 channels), only captures a smoothed representation of the shallow Earth structure. Therefore, V S30 from the velocity model may not fully capture the structural changes that can rapidly occur at shallow depths. Secondly, we expect the ACFs to have different depth sensitivities based on their frequency ranges. Therefore, a single parameter, namely V S30 , does not account for such depth sensitivity variations. Thirdly, while an accurate value of V S30 can be useful for some geotechnical engineering purposes, it may not be the best parameter to explain the intensity of dv/v drops. For example, in a 30 m profile composed of a very shallow and soft sediment layer overlaying a stiffer material, nonlinearity is expected to only occur in the first layer. Therefore, a complete velocity profile of each site is likely to be more informative than a summarizing parameter (i.e., V S30 ; [START_REF] Bonilla | Contribution of k-net and kik-net data to the monitoring of nonlinear properties of the shallow crust[END_REF], and future work should focus on refining the shallow velocity structure.

Earthquake source duration increases as seismic moment rises, which expends the lower frequency content of the recorded seismic waves in the far field [START_REF] Brune | Tectonic stress and the spectra of seismic shear waves from earthquakes[END_REF]. Such frequency content changes can potentially create a time shift of the ACFs first negative arrivals and therefore introduce source-related artifacts in the dv/v measurements. However, earthquake source effects are likely small compared to nonlinear effects in this study for three reasons. First, [START_REF] Qin | Imaging and monitoring temporal changes of shallow seismic velocities at the garner valley near anza, california, following the m7.2 2010 el mayor-cucapah earthquake[END_REF] showed that ACFs calculated from earthquake waveforms recorded at a surface station at a sedimentary site display similar dv/v values as those computed between borehole and surface stations, which remove the contribution of the earthquake source and path effects. This demonstrates that nonlinear effects at sedimentary sites are generally significantly stronger than earthquake source effects. Secondly, larger magnitude earthquakes do not necessarily cause larger ACF time delays (Figure S3). Thirdly, earthquake source effects are minimized by the range of earthquake magnitudes used to compute the reference ACFs. Each reference ACF is obtained by stacking 20 ACFs from earthquake velocity magnitudes primarily ranging between 1.4 and 4.3 (Figure S4), which overlaps the 1.4-5.3 magnitude range of all the earthquakes analyzed in this study. While this processing step likely smooths earthquake source effects, all dv/v measurements are made with respect to a reference waveform, and are therefore not absolute. Yet, earthquake source effects could become significant if reference ACFs are exclusively computed from small magnitude earthquakes.

While the dynamic strains analyzed in this study are relatively weak, we observe significant relative velocity changes from the ACFs, which indicate a nonlinear response of marine sediments. Due to the weak levels of shaking, soil nonlinearity only occurs during the passing of seismic waves, and no long-term effects can be detected. Nevertheless, the nonlinearity thresholds of strain levels obtained along the cable are consistent with those from laboratory experiments [START_REF] Pasqualini | Nonequilibrium and nonlinear dynamics in berea and fontainebleau sandstones: Low-strain regime[END_REF][START_REF] Remillieux | Propagation of a finite-amplitude elastic pulse in a bar of berea sandstone: A detailed look at the mechanisms of classical nonlinearity, hysteresis, and nonequilibrium dynamics[END_REF]Ten-Cate et al., 2004), although larger dv/v drops are observed in this study. This difference could be explained by the different materials being analyzed (e.g., sandstone versus marine sediments), with marine sediments having a higher susceptibility to nonlinearity. Our dv/v results are also comparable to those from (Bonilla & Ben-Zion, 2020) who observed dv/v drops of a few percent during small local earthquakes in California. To further validate our approach, we also compute ACFs from earthquake data recorded by the horizontal accelerometer along the axis of the cable from the SOB3 station (Figure S5). The ACFs from the SOB3 station exhibit similar features, with a nonlinearity threshold of the same order as that obtained with the DAS data, which validates our approach.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Geophysical Research Letters

The dynamic peak strains recorded by the DAS channels are in the direction of the cable; however, DAS has different theoretical sensitivities depending on the type of seismic waves and their incidence angles [START_REF] Martin | Introduction to interferometry of fiber-optic strain measurements[END_REF]. For example, DAS records from earthquakes occurring at a 90-degree angle from the direction of the cable are expected to exhibit less energy than events happening along the axis of the cable. Moreover, the relatively long gauge length (i.e., 40 m) used to record the DAS data could create notches in the frequency spectrum between 2 and 40 Hz [START_REF] Dean | The effect of gauge length on axially incident p-waves measured using fibre optic distributed vibration sensing[END_REF]. Nevertheless, the steep subduction zone in the Tohoku region [START_REF] Hayes | Slab2, a comprehensive subduction zone geometry model[END_REF] combined with shallow and slow sediment layers interfere with the propagation of seismic waves, which likely arrive with almost vertical angles to the cable. This translates into high apparent velocities of all earthquake wavefields recorded by the cable (e.g., Figure 1c), which limits both the azimuthal and gauge length effects on the recorded data. We further confirm this point by comparing the maximum amplitudes of DAS and SOB3 data during the 105 earthquakes considered in this study in Figure S6. Both datasets exhibit similar maximum amplitudes with respect to azimuth angles to the earthquake epicenters, which confirms that there is no noticeable azimuthal effect for the DAS data.

The largest ground motions during the two-week experiment occurred during a M V 5.6, which occurred 50 km east of the SOB3 station (Figure S7). However, the amplitude of the DAS data clipped and are therefore not usable in our analysis. The data clipping issue is caused by rapid phase changes that occurred during strong ground motions, which wrap the signal's phase. To reduce clipping effects and improve the dynamic range

of DAS experiments, one can increase the laser pulse rate frequency, which would limit the maximum distance that the DAS interrogator can sample, and/or reduce the gauge length, which could result in a decrease of the SNR of the recorded wavefield [START_REF] Mellors | An examination of das as a possible earthquake early warning tool[END_REF]. Despite these drawbacks, a better tuning of the DAS parameters could allow us to record strong ground motions that are likely to trigger stronger nonlinear soil responses.

Conclusions

We analyzed the ground motions of 105 earthquakes recorded along a fiber-optic cable during a two-week DAS campaign offshore Japan. We computed ACFs of earthquake ground motions and detected relative velocity changes in the marine sediments surrounding the cable from the ACFs. Large drops of dv/v are observed along the cable and are typical of nonlinear behavior of the medium. Moreover, the dv/v changes are frequency and spatially dependent, highlighting the sensitivity of DAS ACFs to the shallow Earth structure.

This study demonstrates that earthquakes recorded by DAS can be used to characterize soil nonlinearity during ground motions. This could be critical for fiber-optic cables being considered for earthquake early warning purposes, as soil nonlinearity impacts the amplitude and frequency content of the recorded wavefield and could bias rapid magnitude estimations. Finally, the ACF approach could easily be applied to other DAS datasets recorded in populated regions located on top of sedimentary basins, such as Mexico City and Los Angeles, to better characterize seismic hazard.

Figure 1 .

 1 Figure 1. (a) Topographic map of Japan and its surroundings including 105 earthquakes. The red rectangle denotes the region shown in (b). (b) Bathymetric map offshore the Sanriku coast including the seafloor cable observation system. The orange line denotes the buried section of the cable between channels 500 and 9000 and the location of channel 5000 is indicated by the red cross. The white circles and purple inverted triangles show the positions of the accelerometers and tsunami-meters, respectively. The location of the MV 2.5 event (red circle), for which the waveforms are shown in (c), and that of other nearby earthquakes (gray circles) are highlighted. The magnitude scale is the same as in (a). (c) Strain waveforms of the MV 2.5 event bandpass filtered between 2 and 30 Hz between channels 500 and 9000. Waveform amplitudes are clipped for visibility. (d) Strain-rate waveforms of the MV 2.5 earthquake bandpass filtered in different frequency bands at channel 5000. The gray area denotes the time period over which the ACFs are calculated. (e) Amplitude normalized ACFs computed from the waveforms shown in (d).

Figure 2 .

 2 Figure 2. ACFs computed from the 105 earthquakes bandpass filtered between 10 and 20 Hz at channels (a) 5000 and (b) 7000. The amplitude of the data is clipped for visibility. (c-d) Same as (a-b) in the 15-30 Hz frequency band. (e-f) Dynamic peak strains after bandpass filtering the earthquake strain waveforms between 10-20 Hz and 15-30 Hz at channels 5000 and 7000, respectively. The ACFs in (a-d) are sorted by the increasing dynamic peak strain values shown in (e-f).

Figure 4 .

 4 Figure 4. (a) Mean dv/v values computed between weak and strong ground motion ACFs at each channel and each frequency band. (b) Average maximum dynamic peak strain over the 10 largest peak strains at each channel and each frequency band. (c) Relative velocity changes as a function of the filtered dynamic peak strain in the 3-6 Hz frequency band. The error bars represent the one standard deviation from the mean of the dv/v measurements (Figure S2). (de) Same as (c) for the 15-30 and 20-40 Hz frequency bands. (f) Relative velocity changes as a function of the average S-wave velocity within the first 30 m of the ground (VS30) for the 3-6 Hz frequency band. (g-h) Same as (f) for the 15-30 and 20-40 Hz frequency bands. In (c-h), the color-scale corresponds to the channel number.

  Distributed Acoustic Sensing Array at Garner Valley, California. Bull. Seismol. Soc. Am., 107 (2), 603-610. Retrieved from https://doi.org/10.1785/ 0120160168 doi: 10.1785/0120160168

Acknowledgments

We thank Takeshi Akuhara for providing useful information about DAS measurements.

We thank Fujitsu for cooperating with the Earthquake Research Institute (ERI), The University of Tokyo, for the DAS measurement campaigns. We are grateful to the Editor Daoyuan Sun and two anonymous reviewers for providing helpful comments that improved this study. Funding: This project was partly supported by the discretionary budget of the director of ERI. The observations were carried out as part of the Earthquake and Volcano Hazards Observation and Research Program by the Ministry of Education,

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved. Open Research: The data and codes developed to perform the technical analysis and to reproduce most figures are available at https://doi.org/10.5281/zenodo.7044595. [START_REF] Hunter | Matplotlib: A 2d graphics environment[END_REF] and some of the data processing steps have been performed using ObsPy [START_REF] Beyreuther | ObsPy: A Python Toolbox for Seismology[END_REF].

All the Figures are plotted with Matplotlib

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Geophysical Research Letters Sawazaki, K., Sato, H., Nakahara, H., & Nishimura, T. (2006). Temporal change in site response caused by earthquake strong motion as revealed from coda spectral ratio measurement. Geophys. Res. Lett.,33 (21). Retrieved from https://