Ambient‐Noise Wave‐Equation Tomography of the Alps and Ligurian‐Provence Basin

Taking benefit of the AlpArray temporary network and permanent networks in W‐Europe, we construct a 3‐D onshore‐offshore velocity model of the crust and upper mantle using ambient‐noise wave‐equation tomography. We use a frequency‐dependent phase traveltime misfit function in an iterative procedure to refine a recent 3‐D Vs model computed from a Bayesian two‐step ambient noise tomography. Observed waveforms consist of vertical‐component noise correlations from 600 broadband stations in the Alps and surroundings, including ocean‐bottom seismometers in the Ligurian sea. We perform 3‐D inversion in the 5–85 s period range. In the long‐period band (20–85 s), an elastic approximation is considered, while in the 5–20 s band, we account for the effect of water layer in the Ligurian sea by applying a fluid‐solid coupling for acoustic‐elastic waveform simulations. The resulting Vs model enhances the shape and contrast of velocity structures, accounting for 3‐D and finite‐frequency effects. It emphasizes the deep sediments of the Ligurian‐Provence basin and focuses on the low‐velocity anomalies of the crust in the W‐Alps. We obtain a high‐resolution Moho depth map covering the Alps and Ligurian sea. In the W‐Alps, this map confirms the deepening of the European crust following the subduction beneath Adria and the existence of major structures such as the Moho jump beneath the external crystalline massifs and shallow depths associated with the Ivrea Body. It provides further constraints on the deep structure beneath the Ligurian‐Provence basin, regarding the lateral and along‐strike crustal‐thickness variations from the oceanic domain to the conjugate margins.

10.1029/2023JB026776 2 of 28 (e.g., Gueguen et al., 1998;Jolivet et al., 2020).This extension started along the Provence region leading to the Ligurian-Provence basin, and has further spread from east to west and south, resulting in the opening of the Algerian basin, and later, of the Tyrrhenian basin (e.g., Rollet et al., 2002;Séranne, 1999).The opening of the Ligurian-Provence basin was initiated at 35 Ma by a rifting phase between Europe and the Corsica-Sardinia block.The progressive south-eastward roll-back and retreat of the Adria slab below the Corsica-Sardinia domain led to Oligocene stretching of the continental crust followed by continental break-up during the early Miocene, and to the genesis of an oceanic crust between 15 and 11 Ma (e.g., Jolivet & Faccenna, 2000;Séranne, 1999).As a result, three principal domains, mainly identified from active seismic data, describe the present-day geological setting of the Ligurian-Provence basin (Figure 1): (a) two thinned conjugate continental passive margins corresponding to the Ligurian-Provence margin to the north and Corsican margin to the south; (b) an oceanic domain in between, described as an "atypical oceanic crust" by Rollet et al. (2002) because it is thinner than normal oceanic crust; and (c) two transitional domains separating the margins and the oceanic domain, likely made up of a very thin continental crust overlying a thick rift-related corner of magmatic underplating (e.g., Séranne, 1999).
The complex geodynamic context associated with the interaction between the opposing vergence subductions of Europe and the Apennines has resulted in a three-dimensional (3-D), heterogeneous lithospheric structure that is difficult to image by seismology.Consequently, questions remain regarding the role of the 3-D geometry of European continental subduction in the present-day architecture of Western Alps (e.g., Malusà et al., 2021), and the role of the rapid retreat of Apennines subduction in the present-day petrological-lithological composition of the crust beneath the Ligurian-Provence basin (e.g., Rollet et al., 2002).Thus, accurate 3-D seismic imaging of the crust and upper mantle is critically needed to complement the 2-D seismological observations, especially since the region benefits from a fairly wide coverage of seismological stations, including the newly deployed AlpArray land-sea network (AASN; Hetényi et al., 2018).The abovementioned models have been derived using conventional ANT, which involves two steps: (a) 2-D traveltime inversion for group or phase velocity maps and (b) 1-D depth inversion for V s based on local dispersion curves.This procedure has two limitations that may bias the geological interpretations: (a) the ray theory assumption when computing 2-D Rayleigh-wave velocity maps, which is only valid in the high-frequency case (e.g., Cerveny, 2003;Snieder, 1986), and (b) the local 1-D nature of the depth inversion in the second step, which does not account for the 3-D lateral heterogeneity of the medium, thus limiting the velocity model to be pseudo-3-D by construction.The aim of our study is to overcome these limitations in the context of Western Alps and Ligurian-Provence basin by building a self-consistent onshore-offshore 3-D velocity model using a wave-equation tomography (WET) of ambient noise.This allows us to take into account the effect of 3-D structures and of the water layer in the Ligurian sea on the propagation of surface waves.Indeed, wave-equation-based tomographic methods are an alternative of choice to overcome such assumptions, as they naturally accommodate for 3-D heterogeneity and finite-frequency effects, thus providing a more realistic sensitivity kernel for surface waves.These approaches consist of iteratively updating the velocity model by minimizing a misfit function between observed and synthetic waveforms obtained through 3-D numerical modeling of seismic wave propagation.
In the following, we refer to ANT when the two-step method is used to invert the correlations, to FWI when the seismic wave propagation is modeled in 3-D by solving the wave equation and the misfit function involves the waveform that is, phase and amplitude (e.g., Virieux & Operto, 2009, and references herein), and to WET when the misfit function involves only the phase (e.g., Luo & Schuster, 1991;Tape et al., 2010).
With the availability of accurate large-scale seismic wavefield modeling and computing resources, wave-equationbased tomographic methods have been widely applied to the crust and upper mantle (e.g., Beller et al., 2018;Fichtner & Villaseñor, 2015;Tape et al., 2010;Yuan et al., 2014;H. Zhu et al., 2012), often providing higher resolution images than those obtained from ray theory and finite-frequency tomography.However, their application to ambient noise data is challenging.Indeed, noise correlations provide the Green's function of the medium only under certain conditions, for instance, when considering a homogeneous distribution of uncorrelated noise sources (Roux et al., 2005;Wapenaar, 2004), or more generally assuming equipartition of the noise wavefield (Campillo, 2006;Lobkis & Weaver, 2001;Snieder et al., 2007;Weaver & Lobkis, 2001).These conditions are rarely fully met in practice.There are two ways to deal with this difficulty.
The first way consists of taking into account the distribution of the noise sources when computing synthetic noise correlations and their sensitivity kernels.Indeed, the cross-correlations waveforms can be strongly affected by the distribution of noise sources (Fichtner, 2014;Tromp et al., 2010) and by the processing applied to the noise records (Bensen et al., 2007;Fichtner et al., 2016Fichtner et al., , 2020)).In particular, the amplitude of noise correlations depends mostly on the energy and distribution of the noise sources (e.g., Hanasoge, 2014).As a consequence, spatial variations of noise correlation amplitude cannot be interpreted unambiguously as lateral contrasts of attenuation (Stehly & Boué, 2017;Tsai, 2011) but can be used to image the distribution of noise sources (Ermert et al., 2016(Ermert et al., , 2017(Ermert et al., , 2021;;Igel et al., 2021Igel et al., , 2023;;Stehly & Boué, 2017).Therefore, methodological efforts have been made to apply FWI to noise correlations without assuming that they are similar to the Green's function of the medium, but by treating them as self-consistent observables.This requires joint invert for the distribution of noise sources and the Earth structure to take into account the distribution of noise sources when computing synthetic noise correlations and their sensitivity kernels (Fichtner, 2014;Hanasoge, 2014;Sager, Boehm, et al., 2018, Sager, Ermert, et al., 2018;Tromp et al., 2010).However, applying this approach to image the crust and upper mantle at a continental scale remains challenging due to the complexity of the dynamics of microseismic noise sources and the heterogeneity of the medium.
An alternative approach is to apply WET (Luo & Schuster, 1991) rather than FWI to noise correlations.Indeed, while the heterogeneous distribution of noise sources and data processing affect the correlations waveforms, traveltime measurements of the correlations surface-wave fundamental mode are less affected (Delaney et al., 2017;Fichtner, 2014;Froment et al., 2010;Tsai, 2009Tsai, , 2011)).Consequently, in crustal environment, the traveltime residuals depend more on the heterogeneity of the medium than on the source distribution (Froment et al., 2010;Kimman & Trampert, 2010;Yang & Ritzwoller, 2008).It is thus possible to assume that at least a subset of the correlations provide accurately the Green's function surface-wave dispersion.Hence to apply WET, noise correlations are modeled assuming that they are similar to the Green's function, that is, without taking into account the distribution of noise sources.In this case, the traveltime (phase) of surface-wave fundamental mode is inverted rather than the waveform.However, this approach has two main limitations: (a) it does not account for the noise source distribution which can bias the apparent arrival time of surface waves (Delaney et al., 2017;Froment et al., 2010;Tsai, 2009).This is mitigated in practice by selecting only station pairs for which the arrival time of surface waves is similar in the positive and negative correlation time; (b) noise correlations that contain mostly the fundamental mode of surface waves are compared with synthetic correlations that contain all possible modes of propagation.This may affect misfit measurements if the different modes are not well separated in time.Nevertheless, WET applications to noise correlations have recently demonstrated its potential to provide accurate velocity models of the crust and upper mantle that better explain the observations than conventional ANT (Chen et al., 2014;Liu et al., 2017;Lu et al., 2020).
Using the ANT model from Lu et al. (2018) as an initial model, Lu et al. (2020) were the first to apply WET to the alpine region using noise correlations.They demonstrated that WET is relevant to obtain models that better explain the data, and to improve the resolution of ANT models in the Alps.However, their WET model does not cover the Ligurian sea since the offshore part of the AlpArray network was not yet available.We refer the reader to Nouibat, Stehly, Paul, Schwartz, Bodin, et al. (2022) and Nouibat, Stehly, Paul, Schwartz, Rolland, et al. (2022) for more details about the comparison of recent ANT methodologies and models in the region.
In this study, we improve the WET methodology used by Lu et al. (2020).Theses improvements consist of: (a) performing an acoustic-elastic coupled WET to better represent the wave-physics and consider more realistic constraints on velocity structure in the oceanic domain of the Ligurian basin and along the continental margins; (b) inverting seismic data in a broader band (5-85 s) to constrain the shallow part of the crust as well as the upper mantle, employing a hierarchical inversion strategy that avoids possible cycle-skipping issues; and (c) using a random sub-sampling scheme over 185 virtual sources rather than a fixed number of virtual sources.Similar to Lu et al. (2020), we minimize the frequency-dependent phase traveltime differences of Rayleigh waves and tackle the inversion using the SEISCOPE SEM46 code originally developed for exploration scales (Trinh et al., 2019).Taking advantage of the densest seismological coverage in Western Europe including the entire AlpArray network, we perform a large-scale WET to refine the land-sea model by Nouibat, Stehly, Paul, Schwartz, Bodin, et al. (2022) obtained from conventional ANT.
The availability of Rayleigh-wave traveltime measurements across and within the Ligurian sea, and the latest developments of full waveform modeling and inversion in the fluid-solid coupled media (Cao et al., 2022), motivated us to go further in the improvement of WET by considering the effect of the water layer on the 3-D propagation of surface waves from ambient noise correlations.
In an oceanic domain, two different types of surface waves can be recorded on the vertical component.At large periods, when the wavelength of the surface waves is large compared to the thickness of the water layer, the effect of the water on wave propagation is negligible, and the surface waves propagate as Rayleigh waves.At shorter periods, however, the effect of the water layer can no longer be neglected and the surface wave propagates as a Rayleigh-Scholte wave, which is a fluid-solid interface wave.Nouibat, Stehly, Paul, Schwartz, Rolland, et al. (2022) have shown that a water layer thicker than 0.5 km can have a significant effect on the Rayleigh-Scholte wave dispersion at period ≤20 s and should therefore be taken into account in tomographic studies.Thus, several studies have taken the water column into account when inverting the local dispersion curve of Ralyeigh-Scholte waves extracted from noise correlation into a local 1-D V s model (e.g., Carvalho et al., 2022;Guerin et al., 2020;Mordret et al., 2014;Nouibat, Stehly, Paul, Schwartz, Rolland, et al., 2022;Wolf et al., 2021).However, this assumes a strictly 1-D coupling in a laterally homogeneous medium.
In our study, the 3-D effects of the presence of water layer are considered on the surface-wave propagation.For this purpose, we apply the technique of acoustic-elastic coupling in the waveform modeling to deal with the seabed relief.The objective is twofold.First, to investigate the influence of the water layer on the propagation of surface waves between onshore and sea-bottom stations, more specifically Rayleigh-Scholte waves.To our knowledge, this is the first study documenting the influence of the water layer on the 3-D propagation of ambient noise surface waves.Second, to incorporate the fluid-solid coupling in the inversion framework, which makes this study the first ambient noise WET to take into account the effect of the water layer.
The structure of the manuscript is organized as follows.Section 2 describes the data set of ambient noise correlations and the initial velocity model used in the WET.The overall methodology is presented in Sections 3 and 4. Section 3 introduces the iterative process of the WET workflow, while Section 4 is dedicated to 3-D modeling of surface waves, with emphasis on the acoustic-elastic case.In light of the specific case of the Ligurian-Provence basin, we document the importance of the fluid-solid coupling for marine crustal imaging based on ambient noise data.Section 5 presents tomography results and related discussions.In Section 6, we assess the robustness of the resulting 3-D velocity model.

Ambient Noise Data
We use a subset of the data used by Nouibat, Stehly, Paul, Schwartz, Bodin, et al. (2022) and Nouibat, Stehly, Paul, Schwartz, Rolland, et al. (2022) consisting of vertical component cross-correlations computed between 600 broadband stations from all available temporary and permanent networks in the Alpine region (Figures 1, during the period 2015-2019.The station array includes the entire AlpArray Seismic Network (AASN) and the Cifalps-2 experiment.
The offshore component of the data set consists of high-quality correlations between the AASN ocean-bottom seismometers (OBSs) obtained from an iterative procedure involving onshore stations of the AASN and permanent networks.Beforehand, the OBSs daily noise records were cleaned from transients (e.g., glitches) and seafloor-induced noises (compliance and tilt).For more details about this specific processing, readers can refer to Nouibat, Stehly, Paul, Schwartz, Rolland, et al. (2022).
We keep only reliable cross-correlations by applying the following criteria: (a) signal-to-noise ratio (SNR) greater than 5, (b) difference in group velocity measured on the positive and negative parts less than 0.2 km/s, and (c) inter-station distance larger than one wavelength for the maximum period considered in each period band.The SNR is defined as the ratio of the maximum amplitude of the surface wave to the standard deviation of the signal starting after the surface wave window.The final data set includes ∼22 × 10 3 to 55 × 10 3 high-quality cross-correlations depending on the period band (Table 1).For the WET, all stations in Figure 1 serve as receivers, out of which 185 stations are used as virtual sources (white circles).

Initial Model
The initial model consists of the 3-D onshore-offshore V s model by Nouibat, Stehly, Paul, Schwartz, Bodin, et al. (2022) Nouibat, Stehly, Paul, Schwartz, Bodin, et al. (2022) and Nouibat, Stehly, Paul, Schwartz, Rolland, et al. (2022) was obtained from ANT using all available permanent and temporary seismic networks in Western Europe from 2015 to 2019.This model is derived from a hybrid data-driven tomography involving a 2-D transdimensional Bayesian inversion for group-velocity maps and their uncertainties, and a 1-D probabilistic inversion for V s at depth.The onshore part of this model has been validated in the western Alps by comparison with other available geophysical studies, for example, receiver functions (RFs), controlled-source seismic (CSS) experiments, Bouguer anomaly, along the Cifalps-1 and Cifalps-2 profiles (Nouibat, Stehly, Paul, Schwartz, Bodin, et al., 2022;Paul et al., 2022).The offshore part has been validated by comparison with a high-resolution V p cross-section derived from refraction and wide-angle seismic profiling along the axis of the Ligurian basin (Nouibat, Stehly, Paul, Schwartz, Rolland, et al., 2022).However, it remains limited by the assumptions made in the two steps of the inversion, in particular the high-frequency assumption for Rayleigh-wave propagation in the 2-D inversion, and the 1-D assumption in the inversion for V s .Because of the 1-D assumption, they assumed a 1-D fluid-solid coupling in the oceanic domain, which is not suitable for dispersion of Rayleigh-Scholte waves or conversion of Rayleigh waves to Rayleigh-Scholte waves across the margins.The objective of the present study is to overcome these limitations and improve the resolution of the V s model by performing 3-D WET.

Iterative Inversion Scheme
Similar to Lu et al. (2020), we formulate the inverse problem as the minimization of the misfit function (ξ) defined by the differences of frequency-dependent phase traveltime between observed and synthetic vertical-component waveforms of ambient-noise Rayleigh waves where ΔT i (ω, m) is the phase traveltime residual for the ith station pair at frequency ω, given the model parameters (m).We measure ΔT using a multi-taper technique (Tape et al., 2010).We tackle the minimization of (ξ) using a quasi-Newton local optimization method, with the following iterative scheme (Métivier & Brossier, 2016): where α k is the step length estimated by line search strategy and δm k = −B k ∇ξ(m k ) is the model update in which B k is an approximation of the inverse of the Hessian matrix (i.e., the second-order derivative of ξ with respect to model parameters).We use the limited memory version of quasi-Newton methods (l-BFGS; Nocedal, 1980) to estimate B k using gradients at few previous iterations.We update the velocity model by inverting from low to high frequencies, in four narrow period bands: 5-10, 10-20, 20-40, and 40-85 s.Given the dispersive character and depth sensitivity of the Rayleigh-wave phase velocity at these periods, such a hierarchical approach is useful to prevent cycle-skipping issue.This approach first constrains large-scale anomalies in the upper mantle part of the model using the long-period content, and then it constrains small-scale anomalies in its crustal part using the short-period content.
We reconstruct simultaneously V s and V p to account for the sensitivity of the Rayleigh-wave to V p , which becomes non-negligible at short periods (e.g., Eddy & Ekström, 2014;Qiao et al., 2018).Therefore, the relative perturbation of the misfit function (ξ) to perturbations of the two parameters is given by the linear relation:  () =   () ln +   () ln , where    and    are the sensitivity kernels for V s and V p respectively (i.e., misfit gradients with respect to V s and V p ). Due to the computational cost of modeling for multiple sources, we employ a randomized sub-sampling strategy of the ensemble of virtual sources over iterations to achieve optimal path coverage without compromising the quality of the inversion results, as usually applied for exploration scales (Warner et al., 2013).Following the workflow shown in Figure 2, our wave-equation-based tomography is summarized in the main steps bellow.
1. Perform the forward modeling for one subset of 64 virtual sources.2. Measure phase traveltime misfits between observed and simulated waveforms and compute the adjoint sources.3. Perform the adjoint modeling and extract gradients from cross-correlations of incident and adjoint fields.4. Gradients are summed and smoothed using a non-stationary Laplace smoothing filter (approximated by a second-order Bessel filter) to remove high wavenumber artifacts and constrain the inversion (Trinh

3-D Simulations of Surface-Wave Propagation
We perform forward and adjoint simulations using the time-domain wave equation modeling of the SEM46 spectral-element-based code.We perform simulations for the vertical component of Rayleigh waves by applying a point-force at the location of virtual sources with a filtered Dirac delta as source time function (identical for all virtual sources).Therefore, the synthetic waveform (u syn ) at a receiver position for model (m) is the convolution product between the source time function (s) and the synthetic Green's function (G syn ) for the source-receiver pair (e.g., Lu et al., 2020): To compare with synthetic waveforms, we convolve the source time function (s) with the time derivative of the noise cross-correlation (C) for the source-receiver pair, which is a good approximation of the Green's function (G) of the target medium (e.g., Lobkis & Weaver, 2001;Roux et al., 2005;Snieder, 2004;Wapenaar, 2004;Weaver, 2005): The time sampling of simulation varies from 0.04 to 0.0015 s, with the Courant-Friedrichs-Lewy condition satisfied (Komatitsch & Tromp, 1999).Surface waves are extracted in a time window with min[d/5, t max − 1.5 T] and max[d/2, t max + 1.5 T] as lower and upper limits in the 10-20, 20-40, and 40-85 s, and a time window with min[d/4, t max − 2T] and max[d/1.5, t max + 2T] as lower and upper limits in the 5-10 s, where d is the inter-station distance (in km), t max is the arrival time of the maximum of the envelope of synthetic waveform, and T is the maximum period of the considered band.

Elastic Modeling of Rayleigh Waves
We simulate the incident wavefield in the 20-85 s period band by solving the second-order elastic wave equation: where u is the displacement field, f s is the force vector, σ and ɛ are the second-order stress and strain tensors, respectively, ρ s is the solid density, and C is the fourth-order elastic stiffness tensor.We employ a flexible Cartesian-based hexahedral mesh with element size varying from 15 to 6 km in the horizontal direction and from 10 to 4 km in the vertical direction (see Table 1).The wavefield is interpolated inside each spectral element using Lagrange polynomials of order 4. We parameterize Gauss-Lobatto-Legendre (GLL) points by isotropic S-wave velocity (V s ) starting from the model described in Section 2, and (V p , ρ) computed from V s (Brocher, 2005).Finally, we vertically deform the mesh grid for an accurate representation of the complex topography(bathymetry) on the free surface and to handle the Earth's curvature.
Misfit elementary gradients with respect to the medium parameters are obtained from the zero-lag correlation of adjoint and incident wavefields, based on the adjoint-state method (e.g., Plessix, 2006;Tromp et al., 2005).Adjoint wavefields are obtained by injecting the adjoint sources at receiver locations in the adjoint wave equations.The adjoint source for the multi-taper method is given by the derivative of synthetic waveforms weighted by the frequency-dependent phase difference measurements (Tape et al., 2010).Gradients of (ξ) with respect to the solid parameters (ρ s , C IJ ) are obtained by the zero-lag cross-correlations where  ū is the adjoint wavefield vector associated with u, and  ε is the adjoint of the strain field ɛ for model sets Ω s .

Acoustic-Elastic Modeling of Scholte-Rayleigh Waves
Since we also invert phase dispersion data for S-wave velocities in the Ligurian-Provence basin, we have to consider the influence of the seafloor topography on the propagation of surface waves in the 5-20 s period band.Indeed, Scholte waves can be excited as a result of the fluid-solid interaction (e.g., Zheng et al., 2013;J. Zhu et al., 2004).Scholte waves can also be generated from Rayleigh-Scholte mode conversion at the continent-ocean transitions.To further investigate this influence, we account for the effect of the water layer by considering fluid-solid coupling in the 3-D wave equation through the second-order equation system (Cao et al., 2022): ( where the first equation describes the elastic wave propagation in the solid domain (Ω s ) in terms of displacement vector u (same as Equation 3).The second equation describes the acoustic wave propagation in the fluid domain (Ω f ) in terms of scalar displacement potential φ.P f is the pressure source associated with a force vector f f , ρ s is the fluid density, and κ is the bulk modulus of the fluid.The third system gives the boundary conditions along the fluid-solid interface (Γ fs ), describing the interaction between the two domains.In our implementation, we consider to discretize the topography of onshore-offshore transition zones through staircase functions to preserve the configuration of structured mesh used in SEM46.We consider a minimum vertical size of 200 m for elements in the water column, ensuring proper sampling of waves in water.Misfit gradients in the solid domain are given by expressions (4).In the fluid domain, gradient building is performed following the hybrid approach proposed by Cao et al. (2022).Gradients of (ξ) with respect to the fluid parameters (ρ f , κ) are then obtained by the zero-lag cross-correlations: 10.1029/2023JB026776 9 of 28 The gradients on V p and V s are estimated using elementary gradients (4) and ( 6) based on the chain rule.
Figure 3 shows snapshots of the simulated vertical displacement recorded at the surface, in the 5-10 s period band for an onshore source station located in southern France.We performed elastic (Figures 3a1-3c1) and acoustic-elastic modeling (Figures 3a2-3c2) in the initial model.Snapshots for the simulation in the 10-20 s period band are shown in Figure S1 in Supporting Information S1. Figure 4 shows synthetic waveform fits data at an OBS located in the central Ligurian-Provence basin at a depth of ∼2,700 m (location in Figure 4a).Figures 4b1, 4b2, 4c1, 4c2, 4d1, and 4d2 show simulations in the initial model, and Figures 4b3, 4c3, and 4d3 show simulations in the final model.In the elastic modeling case, the source-OBS pair is assumed to be located on the free surface.
We therefore observe the propagation of a Rayleigh wave everywhere in the oceanic part (Figures 3a1-3c1).In the acoustic-elastic case, we account for the fluid-solid coupling by considering the presence of the water layer in the Ligurian sea.In this case, we observe a mode conversion from Rayleigh to Rayleigh-Scholte wave at the Gulf of Lion margin (Figures 3a2-3c2).The resulting wave packet is characterized by stronger amplitude and slower propagation velocity and includes strong wavefront distortions.However, we notice in both cases a dissipation within the Ligurian-Provence basin and almost the same wavefield in the continental part, as expected.For instance, both simulations evidence a scattering of the wave packet to the north of the Ligurian sea (dashed black frame in Figures 2c1 and 2c2).The discrepancies in amplitude and traveltime in the oceanic part are significant.Indeed, the arrival time of the Rayleigh-wave envelope obtained from elastic modeling in the 5-10 s band is ∼35 s earlier than that of the Rayleigh-Scholte wave obtained from acoustic-elastic modeling (Figures 4b1 and 4b2).Furthermore, the arrival time of the observed surface wave is much closer to the arrival time of the simulated signal with fluid-solid coupling (Figure 4b2).This fit is further improved after inversion as shown in Figure 4b3.Differences between the two modeling results are also visible in the 10-20 s period band with a lower but still significant misfit due to the free-surface boundary condition (Figures 4c1-4c3).As illustrated in Figures 4d1 and 4d2, the differences between elastic and acoustic-elastic simulations are negligible in the 20-40 s period band, showing that the surface wave in these period ranges becomes insensitive to the coupling effect and presence of water.Furthermore, the waveform computed from elastic modeling in the final model fits very well the observed signal (Figure 4d3), indicating that the elastic assumption is realistic-enough at long periods.
Figure 5 shows snapshots extracted at the solid surface (topography and seafloor) and at 10 km depth, in the 5-10 s period band for an OBS source located in the northeastern part of the Ligurian-Provence basin.We observe a significant amplitude attenuation at the depth of the wave propagating in the Ligurian-Provence and Tyrrhenian basins.This is consistent with the trapping of the Rayleigh-Scholte wave in the vicinity of the elastic-acoustic interface (e.g., Nayfeh, 1995).Hence, we detect the long-period dispersion of the elastic wave at 10 km depth, characterized by less distorted wavefront.Finally, Figure 5 shows that wave packets are dissipated in the oceanic part (Ligurian-Provence and Tyrrhenian basins) or through the sharp Provence margin toward the continental part.This suggests that considering the fluid-solid interaction would help to recover the shallow 3-D velocity structures, as it will be illustrated in the next section.

Results and Discussion
As Rayleigh waves from ambient-noise correlations in the period range 5-85 s are mainly sensitive to shear-wave velocities at depth, that is, the most significant variations from the initial model would occur in V s , we analyze here depth slices and Moho maps from the obtained 3-D V s model and discuss their geological implications.

Depth Slices
Figure 6 shows the depth slices of the initial and final V s models at 6 km (Figures 6a1 and 6a2) and 26 km (Figures 6b1 and 6b2) and corresponding relative variations (Figures 6a3 and 6b3).We observe a roughly preserved geometry of velocity structures but strong relative variations.At 6-km depth, the final model exhibits lower velocities (≤2.5 km/s) with strong negative velocity variations (≤−8%) in the Ligurian-Provence basin.These low velocities probably reflect the slow propagation of Rayleigh-Scholte waves, induced by the sedimentary cover.Another new pattern appears in Corsica, where the western Variscan Corsica (see Figure 1) exhibits higher velocities (≥3.5 km/s) than the northeastern Alpine Corsica.Such variations are probably due to Rayleigh-Scholte mode conversions across the Corsican margin, which are taken into account by the 3-D fluid-solid coupling.On the onshore domain, relative variations roughly coincide with geological structures, with a velocity decrease in the western and eastern parts of the subduction wedge and an increase in the forelands.
The 26-km maps (Figures 8b1-8b3) display two striking features: (a) a velocity decrease (down to −8%) in the subduction wedge with a northeastward focusing of the low-velocity anomaly (LVZ) beneath the W-Alps, and (b) a velocity increase (up to 8%) in the crust beneath the Alpine forelands and the Apennines.The transitions between these domains show small variations to the initial model, highlighting the effectiveness of the ANT model in dealing with sharp transitions along the main tectonic boundaries (e.g., the subduction wedge-European foreland boundary, see Figure 1).We will see in Section 6.1 that such transitions exhibit a low average misfit.Information S1.The initial and final V p models remain fairly close due to the limited sensitivity of Rayleigh waves to P-wave velocity (Figure S3 in Supporting Information S1).

Vertical Cross-Sections
We now focus on the western Alpine and Ligurian-Provence regions, and compare the initial and final Versus models along 2-D cross-sections of the crust and uppermost mantle (profiles P1-3 in Figure 1).These profiles traverse complex geological domains that are well known at the surface and have already been investigated at depth using other techniques with independent data.The underlying 3-D structure is therefore particularly interesting for validating the continental and oceanic parts of the model.

Cross-Sections in the Western Alps
The profile P1 (Figure 1) coincides with the Cifalps transect, which was extensively investigated in a number of tomographic studies (Malusà et al., 2021, and references herein).The profile P2 is located close to ECORS-CROP wide-angle reflection transect (Nicolas et al., 1990) and Cifalps-2 RFs transect (Paul et al., 2022).The two profiles cross the W-Alps from the European foreland (subducted lower plate) to the Adriatic foreland (upper plate), and the subduction wedge delimited by the PFT to the west and the IF to the east (Figures 7a1-7b1).
The cross-sections in initial (ANT) and final (WET) V s models as well as their relative velocity changes are shown in Figure 7.The final model displays lower velocities in the European mantle and higher velocities in  (Zhao et al., 2015).Note the remarkable agreement with the V s Moho proxy in (a3).ACC, Adria continental crust; ASB, Adria seismic body; ECC, European continental crust (Nouibat, Stehly, Paul, Schwartz, Bodin, et al., 2022).Labels 1, 2, and 3 indicate velocity structures discussed in the text.
Adria mantle (±2-3%).We note a progressive increasing of velocities beneath the subduction wedge, as highlighted by positive-velocity variations from 40 km depth and 180 km offset.Such variations at large depths indicate a significant contribution to the model update from the long-period component of the dispersion data.The European crust (ECC) is marked by an attenuation (dv/v ≥ +5% in P1) of the lower crust LVZ (label "1") beneath the foreland, and strengthening (dv/v ≤ −3% to 4%) of the LVZ (label "2") beneath PFT trace.These two velocity changes are already visible in Figure 6, indicating a north-eastwards shift of LVZ1.This shift is likely the result of correcting for a bias from the seismic ray deviation that is not accounted for in the ANT model due to the high-frequency approximation.It also reflects the higher resolution of the final model.The resolution tests in Section 6.4 show that such perturbations are fairly reliable.The final model exhibits smaller changes in Adria crust (ACC), which is less structurally complex than ECC, indicating that it is indeed well constrained at depth by the ANT.
Similarly to RFs Moho (dashed yellow line), the 4.2 km/s contour of the final model dips continuously following the subduction of European lithosphere, reaching a maximum depth of ∼80 km.In the north (profile P2), the European Moho reaches only 55 km maximum depth, which is consistent with CSS and RFs data along ECORS-CROP and Cifalps-2.We therefore define this contour as proxy for Moho in W-Alps.Below the subduction wedge (profile P1), the Moho geometry is more stable than that of the ANT which extends to depths of more than 90 km (location "3").Similar changes along P2 underline a decrease in the maximum depth (55 km instead of 70 km).We note strong changes in the V s gradient at Moho depths, particularly along P2 (dv/v ≥ +5%).However, no changes are observed at the Moho jump in this profile, indicating that this structure is well constrained by the ANT.Finally, the crust-sediment interface, which can be roughly approximated by the 2.8 km/s contour, shows a velocity decrease below the Po plain and below the Rhone valley.Such variations along the main interfaces are likely the result of correcting for a bias from the 1-D assumption in the construction of the ANT model.Moreover, their sharpness undoubtedly demonstrates the resolving capacity of our WET.

Cross-Section in the Ligurian-Provence Basin
We now focus on the offshore part of our study area and discuss the P3 profile along the axis of the Ligurian-Provence basin (location in Figure 1).This transect crosses the three main domains involved in the genesis of oceanic crust.The oceanic southwestern part of the transect coincides with the seismic refraction and wide-angle reflection line by Dannowski et al. (2020), while the northeastern part coincides with the marine part of the wide-angle reflection line by Makris et al. (1999).
Figure 8 shows that the final model displays major velocity changes with respect to the ANT model: (a) strong decrease in sediments and crust of the oceanic domain (dv/v ≥ −8%); (b) strong increase in sediments and decrease in crust of the transitional domain; and (c) decrease in sediments and crust of the southwestern margin domain, and an increase to the north-est.Such changes at shallow depths are likely related to the consideration of 3-D fluid-solid coupling.Conversely, the uppermost mantle displays lower velocity changes with an overall decrease in the oceanic domain and an increase in the transitional and continental domains.Similarly to Nouibat, Stehly, Paul, Schwartz, Rolland, et al. (2022), we define the Moho proxy as the 4.1 km/s contour which is more consistent with the V p Moho proxy defined as the 7.3 km/s V p contour by Dannowski et al. (2020) (dashed yellow line).Depth and geometry of the Moho proxy are almost unchanged in the oceanic part.However, we observe a significant increase in velocity reflecting a stronger gradient at this interface, similar to Moho in W-Alps.Finally, we observe a shallower Moho in the final model at the transition to the margin.

New Large-Scale Map of Moho Depth
In a tectonically complex setting as in the Alps and Northwestern Mediterranean, a precise definition of the 3-D Moho structure is fundamental to understand the link between mantle dynamics and the geological setting.
The Moho depth cartography in the Alpine region has been the target of a variety of seismological methods, such as interpolating data from CSS, RFs, and LET (Figure 9a; Diehl et al., 2009;Di Stefano et al., 2009;Lombardi et al., 2008;Piana Agostinetti & Amato, 2009;Spada et al., 2013;Waldhauser et al., 1998), conventional ANT (e.g., Kästle et al., 2018;Lu et al., 2018;Molinari & Morelli, 2011) and more recently from WET (Figure 9b; Lu et al., 2020).Zhao et al. (2015) showed the first evidence of continental subduction beneath the western Alps obtained by analyzing RFs along the Cifalps-1 seismic profile.By using simultaneous observations from RFs across the Cifalps-1 and Cifalps-2 profiles and the V s model by Nouibat, Stehly, Paul, Schwartz, Bodin, et al. (2022), Paul et al. (2022) documented a marked change in the depth of the European Moho from north to south of the western Alps, which occurs over a distance of a few tens of kilometers only.
In this section, we present a new 3-D Moho map at the scale of the Alps and Northwestern Mediterranean (Figure 9c).The Moho depth is estimated from the depth of the iso-velocity surface V s = 4.2 km/s.In the western Alps, our 3-D map of the Moho is overall similar to the one obtained by Lu et al. (2020), with some important discrepancies (Figures 9b and 9c).In the southwestern Alps, while Lu et al. (2020)'s model exhibits ∼40-45 km depths, the European Moho in our model deepens to depth greater than 60 km (see Figure 7), which is more consistent with the continental subduction from RF observations (Zhao et al., 2015).Using the iso-velocity 4.2 km/s as a proxy, we obtain a shallower Moho (20-25 km) associated with the Ivrea Body (white dashed line in Figure 9c).The strong along-strike variations in our Moho depth from north to south are in line with RF observations between the northwestern and southern Alps (Paul et al., 2022), indicating a non-cylindricity of the deep structure of the western Alps.Besides these differences, we observe a similar Moho topography to Lu et al. (2020).In particular, we observe a Moho jump below the external crystalline massifs that in our model appear to be linear and aligned with the Variscan Accident (arrows on Figures 9b  and 9c).

Moho Topography Beneath the Ligurian Sea
3-D crustal imaging of oceanic domains based on ambient-noise surface waves is a challenging topic due to the difficulty of accounting for the effects of the water layer and complex seabed relief in a hybrid medium, as illustrated in this study.One of the main objectives of this study is to produce a homogeneous 3-D map of the Ligurian Moho in line with the physics of surface-wave propagation.Unlike the Alpine region, the 3-D structure of Moho in the Ligurian Sea is poorly constrained by 3-D velocity models.In this region, the few constraints on the Moho depth are mainly provided by CSS data (e.g., Contrucci et al., 2001;Dannowski et al., 2020;Makris et al., 1999;Rollet et al., 2002).This gap is partly due to the poor coverage of seismological stations in this area before the deployment of AlpArray OBSs, making it difficult to achieve 3-D high-resolution imaging with tomography methods such as ANT or LET.
Figure 10 shows the Moho depth map in the Ligurian-Provence basin estimated from the iso-velocity surface V s = 4.1 km/s.The Moho map shows strong lateral variations with depths, ranging from 12 km in the axis of the basin to ∼20-25 km in the conjugate margins of Provence and Corsica.The transition to the Ligurian-Provence and Corsica margins is characterized by a strong gradient in Moho depth while the transition to the Ligurian margin is smoother.In the oceanic domain, the Moho gradually deepens along the axis of the basin, from 12 km in the northern and central parts to ∼14-18 km in the southern part.Red dashed shows the depth contours of a recently published Moho surface from conventional ANT (Magrini et al., 2022).This map shows smoother depth variations than ours within the basin and at transitions to the conjugate margins.Along the basin axis, it shows depths increase from the southern part (15 km) to the northern part (18 km) of the basin.While our Moho has a depth of about 12 km in the oceanic domain traversed by the active seismic profile of Dannowski et al. (2020), the Magrini et al. ( 2022)'s Moho is 15-18 km deep in the northern portion of the profile and more than 18 km deep in the southern portion.However, Dannowski et al. (2020) detected a Moho of ∼12 km depth, which is more consistent with the depths displayed by our model (see Figure 8).These comparisons demonstrate that our Moho model provides a better resolution of the crustal thinning beneath the Ligurian-Provence basin.

Model Robustness
In this section, we will document the robustness of our resulting velocity model.First, we analyze the distribution of phase traveltime misfits of surface waves from the ambient-noise database used in this study.Then, we validate our model with earthquake seismic waves in the western Alps, (a) by comparing P-wave, Rayleigh-wave observed waveforms with synthetic waveforms computed in our model, and (b) by analyzing the distribution of first-arrival traveltime residuals between picked P, S arrivals and synthetic arrivals predicted in our model.

Misfit Analysis
We first document the model robustness by analyzing histograms of misfit for the initial and final model and their spatial and azimuthal distributions.The misfit is expressed as traveltime delay for waves propagating on 100-km distance (s/100 km) to avoid biases related to long inter-station distances.Figure 11 shows histograms of traveltime misfits for the initial (in red) and final (in blue) models at 8, 20, 35, and 35 s period.The mean and standard deviation of the misfit distribution allow a quantitative comparison.The initial histograms are rather narrow and globally centered around 0 misfit value (absolute misfit ≤0.1 s/100 km).The initial standard deviation increases with period (0.3 s/100 km at 8-0.73 s/100 km at 55 s).We do not observe a significant shift of the mean value, which illustrates the consistency of our ANT model.Since our initial model does not suffer from biases toward wrong mean-value velocities, our inversion further contributes in refining the intrinsic shape and magnitude of velocity structures rather than the overall mean velocity at depth.The final model fits the data significantly better than the initial model, with overall smaller average misfits and standard deviations.Histograms of cross-correlation type traveltime differences are shown in Figure S4 in Supporting Information S1.
Figure 12 shows misfit maps at 8, 20, 35, and 55 s periods, obtained by averaging misfit values for the initial and final models over 0.2° × 0.2° cells assuming great-circle ray paths.The spatial distribution of the misfit allows  12a2-12d2).Although the strong misfits in the south of the study region are significantly smaller, some residuals remain, in particular along boundaries of the study region where path coverage is poorer (e.g., southern Corsica, in the 20-s and 35-s maps, or in the western Ligurian-Provence basin at 55 s).The azimuthal distribution of misfit for the initial and final models is shown in Text S2 and Figure S5 in Supporting Information S1.

Validation With Earthquake Waveform Modeling
We validate our final (V p , V s ) model by performing 3-D wave propagation simulations of earthquakes in Western Alps.The synthetic waveforms are simulated up to 1 Hz using the same spectral-element-based solver described in Section 4. Since the absolute waveform amplitude can be affected by the attenuation quality factors (Q p , Q s ) at such high frequencies, and we are interested in comparing both amplitudes and phases between observed and simulated waveforms, we consider here a visco-elastic wave propagation.The model is discretized with an adaptive mesh built in the Cartesian system.We parametrize the GLL points with isotropic V p , V s from our model, and isotropic ρ, Q p , Q s computed from these velocities by empirical formulas (Brocher, 2005).We represent the earthquakes as a moment source, using moment tensor solutions obtained by Petersen et al. (2021) from centroid moment tensor (CMT) inversion using the AlpArray network.We perform simulations using a Gaussian source-time function with half duration determined by the seismic moment (Komatitsch & Tromp, 2002).The seismograms are filtered in three narrow period bands (2-5, 3-7, and 5-10 s) and the time windows are defined as [t max − T, t max + T ] for 10.1029/2023JB026776 20 of 28 the P-wave-packet and [t max − 1.5 T, t max + 1.5 T] for the Rayleigh wave-packet, where T is the maximum period of measurement and t max is the maximum envelope arrival time.t max is detected in velocity range for P-waves (5.0-8.0 km/s) and for Rayleigh waves (1.5-3.5 km/s).Finally, we measure the cross-correlation-based time shifts to quantify the mismatch between the synthetic waveforms and observations.Figure 13 shows an example of a comparison between synthetic and observed waveforms for an earthquake of magnitude M w = 4 and 5-km depth that occurred in North Western Alps in 2017 (Figure 13a for location).The seismograms are recorded by six stations from permanent networks (FR, France; CH, Switzerland; and GU, Italy) and the temporary AlpArray network (Z3).We observe a striking fit between synthetic and observed waveforms for both P and Rayleigh waves, as supported by an overall time delay lower than 0.65 s (Figures 13b-13d).As we expect, the Rayleigh waves are quite well fitted in the 5-10 s period band (Figure 13d), which is the shorter period band used in our inversion.Besides that, the striking fit in the 3-7 s period (Figure 13c) highlights the potential of our model to recover the shallower and lateral small-scale structures.The good coherency in phase and amplitude for the P-waves in the 2-5 and 3-7 s period bands (Figures 13b and 13c) highlights the potential of our model to recover the small-scale intra-crustal structures.This example illustrates that our WET model could be a relevant starting model for WET of regional earthquake body and surface waves-a natural extension of the WET of ambient noise correlations.A more comprehensive validation of this model will be the subject of a future publication.

Validation With First-Arrival Traveltimes of Local Earthquakes
We further document the robustness of our final velocity model by comparing observed and synthetic body-wave traveltimes associated with more than 75 × 10 3 earthquakes in the region of Western Alps, over a period of more than 30 years (Figure 12a).We use ∼1.6 × 10 6 picked times (93 × 10 4 P-wave picks and 66 × 10 4 S-wave picks) coming from the data set of Potin (2016) for the period between 1989 and 2014, complemented by picked times available from the European plate observing system (EPOS) for the period 2015-2021.
The velocity models are discretized according to a regular Cartesian grid of 5 km cell size in horizontal directions and 2 km in the vertical direction.Hypocenters are relocated in our velocity models using the NLLOC software (Lomax et al., 2000).Initial first-arrival traveltimes are predicted in our 3-D V s and V p final models through an efficient eikonal solver (Podvin & Lecomte, 1991) on a finer interpolated square grid of 2 km.Using receivers as eikonal sources saves computer resources, and, therefore, back-tracing rays between each couple (source and receiver) are obtained.Finally, traveltimes are estimated along these rays with a ray sampling of 0.5 km in order to achieve higher-computational accuracy.Using the same procedure, we calculate traveltimes in a 3-D model obtained from a local earthquake tomography using the same data set (Virieux et al., 2023).
The P and S residuals from our model (blue histograms in Figures 14b and 14c) show normal distribution shapes that are rather narrow (standard deviation < 1 s/100 km) and remarkably well centered around 0 misfit (mean < 0.15 s/100 km).The standard deviation of the S residuals is slightly larger, but it should be noted that the total number of P picks is almost 1.4 times larger.However, the mean value of the S distribution is smaller.Such small mean values indicate that our model does not suffer from biases toward wrong mean-value velocities.These histograms obtained from a completely independent data set are comparable to those in Figure 11 based on surface waves from ambient noise.To further assess the relevance of our model, we compare with residuals from the LET model (red histograms in Figures 14b and 14c).The distribution of residuals from our model is notably narrower than that of the LET model.The overall good agreement between observed first-arrival traveltimes and those predicted in our model indicates that it is an accurate enough initial model for 3-D local earthquake tomography.Since it has been obtained using ambient-noise surface waves, this model has the advantage to constrain the lower crust that is often poorly sampled by earthquakes.

Resolution Tests
We perform 3-D spike tests to evaluate the resolution of the final V s model.Three spikes of ∼25 km radius and ∼20 km thickness are located in the crust, beneath Western Alps, Central Alps, and northern Ligurian-Provence basin (Figure 15a1).A spike of ∼50 km radius and ∼40 km thickness is located in the uppermost mantle beneath the northwestern Po basin (Figure 15a2).These patterns correspond to velocity perturbations of up to ∼8% in the crust and ∼5% in the uppermost mantle, with respect to the initial model.Using the same station-pair coverage as our observations, we apply the methodology described in Section 3. Inversion results are shown in the right-hand panel of Figure 15.
At 20 km depth, the crustal anomalies are well recovered.Overall, the shape of the anomalies appears to be slightly stretched.At 50-km depth, the perturbation is recovered but its shape is affected by a strong eastward horizontal smearing.The peak amplitude is weaker as for the crustal anomalies.In both cases, the recovered peak is approximately in the right position.The cross-sections (Figures 15a3-15b3) show that most of the smearing occurs in the horizontal direction.The discrepancies between input and recovered patterns are mainly due to path coverage, to the propagation of Rayleigh wave which induces lateral smoothing of heterogeneities in the crust and upper mantle, and to its sensitivity to the medium which decreases with depth.Nevertheless, part of these discrepancies is also attributable to the tomographic method itself and the convergence rate.
Though these spike tests do not stand for the entire study area, they indicate that our model is fairly resolved in the crust and can be used for the geological interpretation of the 3-D geometry of crustal structures and for their petrophysical characterization.While the resolution in upper mantle remains tolerable for a first-order geological interpretation, it is obviously not sufficient for precise petrophysical characterization.

Conclusion
Using ambient noise data from exhaustive coverage of permanent and temporary arrays in the western European region, we derive a 3-D onshore-offshore velocity model covering the Alps and Ligurian-Provence basin, using WET.We iteratively refine the recent ANT model of Nouibat, Stehly, Paul, Schwartz, Bodin, et al. (2022) by minimizing the phase traveltime differences between observed and simulated Rayleigh-wave waveforms in the 5-85 s period band.Observed signals are obtained from ambient noise cross-correlations and synthetics are computed from SEM-based 3-D elastic and acoustic-elastic modeling of surface waves.The specificity of this study is to highlight the effect of the water layer on the 3-D propagation of Rayleigh waves by applying a fluid-solid coupling for 3-D acoustic-elastic simulations, taking the Ligurian sea domain as an illustration.We demonstrate that the elastic propagation assumption is no longer valid at short periods (5-20 s), since the surface-wave packet is dominated by a composite Rayleigh-Scholte mode propagating with lower velocities.Finally, we incorporate the fluid-solid coupling in the inversion framework of the WET.
In line with the true physics of surface-wave propagation, the WET corrects for the biases in the ANT model related to the high-frequency assumption and the 1-D inversion.The resulting model has a better resolution, with significant intra-crustal changes.In the superficial part, the WET better emphasizes the sediments of the Ligurian-Provence domain and a high-velocity anomaly beneath the Variscan part of Corsica.This improvement is partly due to accounting for Rayleigh-Scholte wave 3-D sensitivity kernels at short periods.In the crust, the WET mainly changes the velocity contrast with an overall velocity decrease in the subduction wedge, focusing the LVZ in the W-Alps to the northeast, and an overall increase in the two forelands.In addition, the WET validates major structures already present in the initial model, such as the subduction of the European crust beneath Adria in the southwestern Alps, and the slow anomalies at the base of the crust in the northwestern Alps, thus validating recent first-order interpretations on the deep structure in this area (e.g., Nouibat, Stehly, Paul, Schwartz, Bodin, et al., 2022;Paul et al., 2022).Significant changes are observed in the oceanic crust of the Ligurian-Provence domain, where intra-crustal velocities decrease significantly along the axis of the basin.
We present a new depth map of the Moho proxy of Western Europe.In the western Alps, this model shows the deepening of the European crust and confirms the Moho jump under the outer crystalline massifs (Lu et al., 2020) and the high-velocity anomalies associated with the Ivrea Body (Nouibat, Stehly, Paul, Schwartz, Bodin, et al., 2022).This proxy is the first truly 3-D representation of the land-sea Moho in the Alps and Northwestern Mediterranean.We show a strong deepening of the Moho toward the Ligurian-Provence and Corsica conjugate margins (from 12 to 20-25 km), and a smoother deepening toward the Ligurian margin and the southern part of the basin (from 12 to 14-18 km).

Figure 1 .
Figure 1.Geological and tectonic setting of the study area, modified from Handy et al. (2010) and Rollet et al. (2002), with locations of seismic stations used in this work (white circles: virtual sources; green circles: receivers).CM, Corsican margin; LPB, Ligurian-Provence basin; LPM, Ligurian-Provence margin.Black lines show the locations of the seismic profiles discussed in the text.
, 2017).The filter is defined by coherent lengths that are directly computed from the inverted parameters V p and V s : L z = L x = L y = 0.1 × λ l where λ l is the local wavelength at the maximum frequency considered. 5. Estimate the step length and update the current velocity model.6. Change the subset for each three iterations.The output model after nine iterations is used as a starting model for the next shorter period band following the same iterative procedure.The density model is updated accordingly using the empirical formula fromBrocher (2005):  = 1.6612 − 0.4721 2  + 0.0671 3  − 0.0043 4  + 0.000106 5  .

Figure 3 .
Figure 3.Comparison of snapshots of the vertical displacement wavefield in the 5-10 s period band, simulated using (a1-c1) the elastic wave equation, (a2-c2) the acoustic-elastic coupled wave equation, for the source station indicated by the red triangle.The wave fields are extracted on a 3-D surface corresponding to the topography in the onshore part and to the bathymetry of the seafloor in the marine part.Black dashed frames show scattered Rayleigh-wave packets.The thin black lines are geological and tectonic boundaries from Figure 1.

Figure 5 .
Figure 5. Snapshots of the wavefield computed from the acoustic-elastic coupled wave equation in the 5-10 s period band, showing the depth attenuation of the Rayleigh-Scholte wave generated by an offshore source (red triangle).Wavefields are extracted (a1-c1) at the surface (Z topo , same definition as in Figure 3) and (a2-c2) at 10 km depth.Black dashed frames show scattered Rayleigh-wave packets.

Figure 7 .
Figure 7.Comparison of depth sections for the initial and final models along the WSW-ENE transect P1 (left-hand side) and the WNW-ESE transect P2 (right-hand side).Locations of P1 and P2 are indicated in Figure 1.(a1-b1) Topographic profiles and geological maps (extracted from Figure 1).(a2-b2) Shear-wave velocities from the ANT initial model (m ANT ).(a3-b3) Shear-wave velocities from the WET final model (m WET ).(a4-b4) Relative variations to the initial model.The white and black curved arrows highlight the subduction of the European lithosphere beneath Adria.The dashed yellow line in (a3) shows the European Moho picked from a receiver functions section(Zhao et al., 2015).Note the remarkable agreement with the V s Moho proxy in (a3).ACC, Adria continental crust; ASB, Adria seismic body; ECC, European continental crust(Nouibat, Stehly, Paul, Schwartz, Bodin, et al., 2022).Labels 1, 2, and 3 indicate velocity structures discussed in the text.

Figure 8 .
Figure 8.Comparison of depth sections for the initial and final models along the SW-NE offshore transect P3 (see location in Figure 1).(a) Shear-wave velocities from the ANT initial model.(b) Shear-wave velocities from the WET final model.The yellow dashed line represents the Moho proxy from the active seismic P-wave velocity section of Dannowski et al. (2020).(c) Relative variations to the initial model.

Figure 9 .
Figure 9. Depth maps of Moho and Moho proxy in the Alps and Northwestern Mediterranean.(a) Spada et al. (2013)'s map composed of three blocks: European Moho (M1), Adriatic Moho (M2), and Ligurian-Sardinia-Corsica-Tyrrhenian Moho (M3).(b) Lu et al. (2020)'s depth map of the V s = 4.2 km/s iso-velocity surface from their WET model.(c) Depth map of the V s = 4.2 km/s iso-velocity surface from our WET model.For a better illustration, the color scale is saturated at 60 km depth; Our model shows 70-75 km depths in the southwestern Alps (continental subduction; see Figure 7).Ad, Adriatic Moho; CA, Central Alps; EA, Eastern Alps; Eu, European Moho; Li, Ligurian Moho; WA, Western Alps.Arrows in purple: Moho step beneath the external crystalline massifs (Lu et al., 2018).Dashed white line: shallow Moho associated with the Ivrea Body.The gray areas hide regions where Moho maps are unconstrained.

Figure 10 .
Figure 10.Depth map of the V s = 4.1 km/s iso-velocity surface from our WET model as a proxy of Moho in Northwestern Mediterranean.Black lines are iso-depth contours from our Moho proxy.Red lines: iso-depth contours from Magrini et al. (2022).Dashed yellow line shows the trace of Dannowski et al. (2020)'s CSS profile.Dashed white lines are limits of the main geological domains in the Ligurian-Provence basin from Rollet et al. (2002).LM, Ligurian margin; LPM, Ligurian-Provence margin; PM, Provence margin.

Figure 11 .
Figure 11.Comparison of histograms of misfit (in s/100 km inter-station distance) for the initial (m ANT , red) and final (m WET , blue) velocity models at 8, 20, 35, and 55 s periods.Labels "mean" and "sd" refer to the mean value and standard deviation of the misfit distribution, respectively.

Figure 12 .
Figure 12.Comparison of the spatial distributions of misfit for the initial (left-hand side: a1-d1) and final (right-hand side: a2-d2) velocity models at 8, 20, 35, and 55 s periods.Misfits are averaged over cells of 0.2° × 0.2° assuming great circle ray paths.The gray area hides regions where the initial model is unconstrained.

Figure 13 .
Figure 13.Comparisons between simulated vertical-component seismograms (red) in the WET model and observed waveforms (black), filtered in the 2-5, 3-7, and 5-10 s period bands for the earthquake and receivers in (a) (red: earthquake, yellow: receiver).The gray areas indicate the P-wave (P) and Rayleigh-wave (R) windows, and ΔT corresponds to the cross-correlation-based time shift between observed and synthetic windowed waveforms (ΔT P and ΔT R , respectively for P-wave, Rayleighwave time shifts).

Figure 14 .
Figure 14.(a) Location map of the seismic events (red dots) and stations (white triangles) constituting the database of first-arrival traveltimes described in Section 6.3.(b, c) Histograms of the traveltime residuals for the P-waves (b) and the S-waves (c) between observations and predicted synthetic arrivals from our velocity model (m WET , blue) and from the LET model (m LET , red).Labels "mean" and "sd" refer to the mean value and standard deviation of the misfit distribution, respectively.

Table 1
and Nouibat, Stehly, Paul, Schwartz, Rolland, et al. (2022), and a V p model converted from this V s Spectral Elements Number/Size and Number of Measurements Used in Each Period Band