A Nouibat
email: ahmed.nouibat@univ-grenoble-alpes.fr

L Stehly

A Paul

S Schwartz

Y Rolland

T Dumont

W C Crawford

R Brossier

Cifalps Team

Alparray Working Group

Ambient-noise tomography of the Ligurian-Provence Basin using the AlpArray onshore-offshore network: Insights for the oceanic domain structure

Efficient processing scheme to remove transients and reduce tilt and compliance from continuous ambient noise recorded by OBSs

• Computation of iterative correlations between OBSs based on a virtual reconstruction of the Rayleigh waves

• Thin, anomalous oceanic crust with gabbroic intrusions evidenced in the basin axis from a joint interpretation of V s and V p models

Plain Language Summary

The Ligurian-Provence basin (Northwestern Mediterranean Sea) is one of the Miocene-Pliocene back-arc basins that resulted from the retreat of the Adria subduction in the plate reorganization due to Africa-Europe convergence. The crustal structure of the basin is still debated, even though it has been probed by active seismic profiling. We compute a high-resolution shear-wave velocity model of the Ligurian-Provence basin and its margins by making optimal use of ambient-noise recordings of seafloor broadband seismometers. In particular, we improve the usually low quality of surface-wave signals in noise correlations between seafloor stations by involving correlations with land stations. The joint interpretation of our S-wave velocity model with a P -wave velocity section obtained in the basin axis by controlled-source seismic profiling provides compelling evidence for the presence of a thick sediment pile above a thin, ∼ 4.5 km-thick oceanic crust, intruded by gabbroic bodies emplaced at the crust-mantle transition. These results show the potential of a joint 1 Introduction

In the last two decades, the deployment of extensive and dense seismic networks of temporary broadband sensors (e.g., USArray, IberArray) has provided a better understanding of the deep structures of the crust and upper mantle, in particular through the emergence of increasingly more precise 3-D seismic imaging (e.g., [START_REF] Moschetti | Seismic evidence for widespread western-us deep-crustal deformation caused by extension[END_REF][START_REF] Levander | Continuing colorado plateau uplift by delaminationstyle convective lithospheric downwelling[END_REF]. However, many of these dense arrays included exclusively onshore sensors, thus preventing 3-D imaging of continent-ocean transitions and oceanic domains. The AlpArray seismic network (AASN) that covers the European Alps and their foreland, is one of the few dense seismic networks consisting of both onshore and offshore stations (Fig. 1; [START_REF] Hetényi | The AlpArray Seismic Network: A Large-Scale European Experiment to Image the Alpine Orogen[END_REF]. The onshore AASN has been used in ambient-noise imaging studies at the scale of the Alps and Apennines to construct 3-D models of shear-wave velocity [START_REF] Lu | High-resolution surface wave tomography of the European crust and uppermost mantle from ambient seismic noise[END_REF][START_REF] Lu | Imaging Alpine crust using ambient noise wave-equation tomography[END_REF], attenuation (coda-Q, [START_REF] Soergel | Coda-Q in the 2.5-20 s period band from seismic noise: application to the greater Alpine area[END_REF], and radial anisotropy (Alder et al., 2021), as well as at the regional scale, e.g., Vienna basin [START_REF] Schippkus | Ambient-noise tomography of the wider Vienna Basin region[END_REF], Western Alps [START_REF] Malusà | The deep structure of the alps based on the cifalps seismic experiment: A synthesis[END_REF], Southeastern Alps [START_REF] Sadeghi-Bagherabadi | High-resolution crustal s-wave velocity model and moho geometry beneath the southeastern alps: New insights from the swath-d experiment[END_REF] and Bohemian Massif [START_REF] Kvapil | Transversely isotropic lower crust of variscan central europe imaged by ambient noise tomography of the bohemian massif[END_REF]. The German-French ocean-bottom seismometers (OBS) of the AlpArray network have been deployed to gain insights into the 3-D structure of the lithosphere beneath the Ligurian-Provence basin and its margins, where major geological/geodynamical issues remain to be clarified. The OBS recordings have been used together with data of onshore permanent and temporary stations in the transdimensional ambient-noise tomography of [START_REF] Nouibat | Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps[END_REF] that covers a large part of Western Europe. In [START_REF] Nouibat | Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps[END_REF], we described the probabilistic inversion strategy for Rayleigh-wave group-velocity maps and their uncertainties, the injection of these uncertainties in the inversion for the V s model, and the validity of the model in the southwestern Alps. The present paper is a complement to [START_REF] Nouibat | Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps[END_REF] that focuses on the Ligurian-Provence basin. We describe here the specific processing of OBS records and the use of iterative noise correlations that are required to improve ray coverage in the Ligurian Sea, hence improving the model resolution in the basin. We further explain how the water layer and its thickness changes are taken into account in the inversion for S-wave velocity. Finally, we compare our V s model to a V p section derived from controlled-source seismic profiling along the basin axis.

Accepted Article

This article is protected by copyright. All rights reserved. The geodynamic context of the western Mediterranean region is controlled by the northward motion of Africa (1 cm/yr) with respect to Europe since Late Cretaceous times (Fig. 1). This global convergence was accommodated by several collision episodes involving Europe with continental micro-plates (Iberia and Adria), leading to the formation of perimediterranean mountain belts (Alps, Apennines, Pyrenees, Dinarides and Betics;e.g., van Hinsbergen et al., 2020). During Miocene and Pliocene times, part of this convergence was accommodated by development of back-arc extensional basins behind the Adria northwestdipping subduction zone (e.g., [START_REF] Gueguen | On the post-25 ma geodynamic evolution of the western mediterranean[END_REF]Jolivet et al., 2020). This extension started in the Ligurian-Provence basin, and has further spread from west to east, resulting in the opening of the Algerian basin, and later, of the Tyrrhenian basin (e.g., [START_REF] Séranne | The gulf of lion continental margin (nw mediterranean) revisited by ibs: an overview[END_REF][START_REF] Rollet | Back arc extension, tectonic inheritance, and volcanism in the ligurian sea, western mediterranean[END_REF]. Crustal thinning in the central Ligurian basin resulted in the formation of a narrow oceanic domain mainly identified from geophysical data, including seismic reflection, refraction and wide-angle profiling (e.g., [START_REF] Egger | The gross structure of the crust under corsica and sardinia[END_REF][START_REF] Déverchère | Malis cruise, rv le nadir[END_REF][START_REF] Dannowski | Seismic evidence for failed rifting in the ligurian basin, western alpine domain[END_REF], altimetry data and magnetic data from aeromagnetic surveys and reduction to the pole [START_REF] Sandwell | Gravity anomaly profiles from ers-1, topex and geosat altimetry[END_REF]. Although the crust of the central Ligurian basin is considered "atypical" because it is thinner than normal oceanic crust and highlights non-linear magnetic anomalies and a concomitant low gravity anomaly [START_REF] Bayer | Magnetic anomaly pattern in the western mediterranean: Earth and planetary sci[END_REF][START_REF] Sandwell | Marine gravity anomaly from geosat and ers 1 satellite altimetry[END_REF][START_REF] Rollet | Back arc extension, tectonic inheritance, and volcanism in the ligurian sea, western mediterranean[END_REF], its petrological and lithological nature is poorly constrained and still debated. Several hypotheses have been considered: (1) thin oceanic crust with tholeitic volcanism overlying mantle rocks, similar to the Tyrrhenian sea (e.g., [START_REF] Mascle | A revised seismic stratigraphy of the tyrrhenian sea: implications for the basin evolution[END_REF][START_REF] Bonatti | Peridotites drilled from the tyrrhenian sea, odp leg 107[END_REF]; (2) partly serpentinized peridotites of exhumed upper mantle mostly devoid of volcanic crust (e.g., [START_REF] Boillot | Undercrusting by serpentinite beneath rifted margins[END_REF][START_REF] Beslier | Ocean-continent boundary in the iberia abyssal plain from multichannel seismic data[END_REF]Jolivet et al., 2020); (3) thinned and stretched continental crust related to an hyper-extended margin (e.g., [START_REF] Mckenzie | Some remarks on the development of sedimentary basins[END_REF][START_REF] Pascal | The ocean-continent boundary in the gulf of lion from analysis of expanding spread profiles and gravity modelling[END_REF][START_REF] Dannowski | Seismic evidence for failed rifting in the ligurian basin, western alpine domain[END_REF].

We combine the OBS records of the AlpArray network with those of 890 onshore stations from the AlpArray temporary network, the Cifalps-2 temporary experiment and European permanent networks (Fig. 1). Ambient-noise data recorded by OBSs have already been used to build shear-wave velocity and anisotropy models in different regions, e.g., Southcentral Pacific [START_REF] Harmon | Using ambient seismic noise to determine short-period phase velocities and shallow shear velocities in young oceanic lithosphere[END_REF], Southeast of Tahiti Island [START_REF] Takeo | Seismic azimuthal anisotropy in the oceanic lithosphere and asthenosphere from broadband surface wave analysis of obs array records at 60 ma seafloor[END_REF] and Western Indian Ocean [START_REF] Hable | Tomography of crust and lithosphere in the western Indian Ocean from noise cross-correlations of land and ocean bottom seismometers[END_REF]. The originality of our approach lies in the specific processing of OBS records and their use within an innovative tomographic framework based on data-driven Bayesian inversions, which has recently been successfully applied to image the lithosphere at the scale of Western Europe [START_REF] Nouibat | Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps[END_REF]. Specific and

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Journal of Geophysical Research: Solid Earth careful pre-processing of OBS records is compulsory because they are affected by noise sources at the seabed such as compliance and tilting, and they may also be impacted by intrinsic instrumental noise [START_REF] Crawford | Estimating shear velocities in the oceanic crust from compliance measurements by two-dimensional finite difference modeling[END_REF][START_REF] Crawford | Identifying and removing tilt noise from low-frequency (¡ 0.1 hz) seafloor vertical seismic data[END_REF][START_REF] Deen | First Observation of the Earth's Permanent Free Oscillations on Ocean Bottom Seismometers[END_REF]. Furthermore, OBSs are sensitive to local noise sources such as tides and currents, boat traffic, or marine animals, which are not recorded coherently over long distances (e.g., [START_REF] Batsi | Nonseismic signals in the ocean: Indicators of deep sea and seafloor processes on ocean-bottom seismometer data[END_REF]. Such noises are therefore unsuitable for ambient-noise tomography, and they even alter the signal-to-noise ratio of surface waves reconstructed by noise correlation.

Finally, the water column above the seismometers (water depth of 1100 -2800 m in our case) may have a significant impact on the quality of the Rayleigh-wave dispersion measurements between distant OBSs, particularly at periods shorter than 15 s. Wolf et al. ( 2021) have highlighted the difficulties in using the AASN OBSs to measure Rayleigh waves dispersion curves from noise correlations. These difficulties are partly related to the high level of seafloor noises. To overcome these difficulties, we propose an innovative way of computing seismic noise correlations using OBS data that consist (1) in a pre-processing of the OBSs noise records that decrease efficiently the seabed noise, and (2) in computing iteratively noise correlations between OBS stations using onshore stations as virtual sources.

The overall methodology and its results are presented in sections 2 to 4. Section 2 is dedicated to the description of a specific pre-processing that aims at cleaning OBS daily noise records from instrumental transient glitches and seafloor noises. Section 3 presents how iterative noise correlations are computed between OBSs using onshore stations as virtual sources. In section 4, we show how noise correlations are used to build a 3-D S-wave velocity model of the Ligurian-Provence basin by computing 2-D probabilistic transdimensional Rayleigh wave group velocity maps and their uncertainties at different periods, and by inverting local dispersion curves to derive a probabilistic V s model. The fifth and final section is dedicated firstly to a validation of the resulting 3-D V s model through a comparison with a recent V p model obtained along a linear refraction, wide-angle seismic profile in the center of the basin by [START_REF] Dannowski | Seismic evidence for failed rifting in the ligurian basin, western alpine domain[END_REF]. Finally, we show how the combination of the V p and V s models along the same profile provides insightful clues to the structure and nature of the crust in the central Ligurian-Provence basin.

Accepted Article

This article is protected by copyright. All rights reserved. We processed ambient-noise records from 23 ocean-bottom-seismometers (OBS in Fig. 1) that were deployed for eight months in the northwestern Mediterranean Sea as the offshore part of the AlpArray temporary seismic network (AASN; network code Z3). All OBSs were deployed in June 2017 by the AlpArray-OBS cruise of the French R/V "Pourquoi-Pas?" [START_REF] Crawford | AlpArray LEG1 cruise[END_REF], and recovered in February 2018 by the MSM71 cruise of the German R/V "Maria S. Merian" [START_REF] Kopp | Rv maria s. merian fahrtbericht / cruise report msm71 lobster: Ligurian ocean bottom seismology and tectonics research, las palmas (spain)[END_REF].

Sixteen LOBSTER instruments (Long-term Ocean Bottom Seismometer for Tsunami The drift of the sensor clocks was measured before the deployment and after the recovery. Recordings were then corrected for clock drift in the RESIF and GEOFON datacenters that archive and distribute the dataset. [START_REF] Hable | Clock errors in land and ocean bottom seismograms: high-accuracy estimates from multiplecomponent noise cross-correlations[END_REF] has demonstrated that the assumption of a linear clock drift is adequate for the ocean-bottom instruments used here.

Glitch removal

The first step of the processing of OBS data is the removal of instrumental transient nearly-periodic impulsive noises (glitches) from continuous records of the 7 French broadband OBSs (BBOBS in section 2.1). Indeed, the vertical-component signals exhibit glitches of 1-hr period caused by the activation of the hourly-check of the internal mass centering of the sensor. The pressure component exhibits glitches with a period of 2.65 hrs, related to data writing on the hard disk. No glitch was detected on the horizontal components. Similar glitches have been observed in other datasets recorded with the same instruments (e.g., [START_REF] Deen | First Observation of the Earth's Permanent Free Oscillations on Ocean Bottom Seismometers[END_REF].

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Journal of Geophysical Research: Solid Earth Similarly to [START_REF] Deen | First Observation of the Earth's Permanent Free Oscillations on Ocean Bottom Seismometers[END_REF], we remove glitches from the data using an average glitch waveform matching algorithm. Detailed explanation can be found in Text S1 of the supporting information. Supplementary Fig. S1 shows a vertical-component daily record of OBS A416A before and after removing the hourly glitches. It documents the efficiency of the processing by comparing the original signal (in blue), the synthetic glitch signal time series (in red), and the final, glitch-free signal (in black). Similarly, Supplementary Fig. S2 shows a daily record of the pressure component before and after removing the 2.65-hrs-period glitches. The power spectral density curves of an example of 1-day raw and pre-processed vertical component record displayed in Fig. 2a, show that removing the glitches reduces the noise level by up to 15 dB at frequencies lower than 10 -1 Hz.

Seafloor-noise reduction

At frequencies below 5 × 10 -2 Hz, the power spectral densities (PSD) of the threecomponent records of all OBSs are dominated by noise due to compliance and tilt (Fig. 2).

The compliance is a long-period pressure signal generated by infra-gravity waves induced by pressure variations in the water column [START_REF] Crawford | Estimating shear velocities in the oceanic crust from compliance measurements by two-dimensional finite difference modeling[END_REF]. Tilting corresponds to displacements and rotations of the sensor induced by seafloor currents that also generate long-period noise. As documented by Fig. 2a, long-period noise is stronger on the horizontal components than on the vertical-component. [START_REF] Crawford | Identifying and removing tilt noise from low-frequency (¡ 0.1 hz) seafloor vertical seismic data[END_REF] have shown that noise on the vertical-component is lower when the instrument is better levelled. However, long-period noise on the vertical component strongly increases if the instrument is tilted, even slightly, since acceleration induced by seafloor currents on the horizontal components is projected onto the vertical component.

Tilt and compliance noises on the vertical-component records are reduced by using a frequency-dependent response function method [START_REF] Crawford | Estimating shear velocities in the oceanic crust from compliance measurements by two-dimensional finite difference modeling[END_REF][START_REF] Crawford | Identifying and removing tilt noise from low-frequency (¡ 0.1 hz) seafloor vertical seismic data[END_REF]. The horizontal components are firstly corrected for compliance noise by subtracting coherent signals derived from the pressure component. In a second step, tilt noise is reduced in the vertical-component by subtracting coherent signals derived from the compliancecorrected horizontal components (black to green in Fig. 2a). Finally, the resulting verticalcomponent signal is corrected for compliance by subtracting the coherent signal derived from the pressure-component signal (green to purple). Text S2 of the supporting information provides a detailed explanation of this procedure.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Journal of Geophysical Research: Solid Earth Figure 2a shows the power spectral density of a daily, vertical-component record of OBS A416A before any correction is applied (blue curve), after correcting from the hourly glitches (black), from the tilt (green) and from the compliance (purple). As shown in Fig. 2b, the coherence of the raw vertical-component with the horizontal components (black and blue curves) increases below 5 × 10 -2 Hz (see also Supplementary Fig. S3a) while coherence with the pressure component increases around 1.2 × 10 -2 Hz. Figure 2a shows that the reduction of the tilt noise (green curve) is maximum (10 dB) at 5×10 -3 Hz, which is the frequency with the maximum coherence between the vertical component and the 2 horizontal components.

The correction for the compliance noise is almost negligible for this record except around 1.2 × 10 -2 Hz, where the coherence between the vertical and the pressure components increases (green area in Fig. 2b). Supplementary Fig. S3b confirms that the correction for compliance is maximized for this station around 1.2 × 10 -2 Hz. In this specific case, the tilt noise is stronger than the compliance noise, which indicates that pressure variations generate less noise than sea-current induced tilt. This strong tilt noise may result from strong currents at the seabed and/or the presence of poorly consolidated sediments directly under the OBS. The lower compliance noise may be due to the large water depth of OBS A416A (∼ 2630 m), damping the effect of pressure-induced infra-gravity waves.

Computation of noise correlations and group-velocity measurements

Once OBS records have been corrected for glitches and seafloor noises, we apply the same pre-processing scheme described below to the records of all stations, onshore and offshore, in order to prepare the calculation of inter-station cross-correlations. As explained in [START_REF] Soergel | Coda-Q in the 2.5-20 s period band from seismic noise: application to the greater Alpine area[END_REF] and [START_REF] Nouibat | Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps[END_REF], we first down-weight the contribution of earthquakes and other high-amplitude transients by removing all 4-hr segments with a peak amplitude 4 times greater than the standard deviation of the current daily record, and with a RMS greater than 1.5 times the daily mean RMS. Each daily record is then filtered into six period bands (3 -5, 5 -10, 10 -20, 20 -40, 40 -80 and 80 -200 s) and amplitudes are normalized by their envelope. Finally, the 6 filtered and normalized signals are stacked to obtain the 4-hr pre-processed broadband signal.

First-order correlations

As in [START_REF] Nouibat | Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps[END_REF], we compute seismic noise cross-correlations for all station pairs by segments of 4 hrs. The 4-hr correlations are normalized and stacked to obtain a

Accepted Article

This article is protected by copyright. All rights reserved. S1a); At frequencies > 5 × 10 -2 Hz, the PSD is dominated by the primary (∼ -140 dB) and secondary (∼ -120 dB) microseisms. At longer periods, the PSD is dominated by the effect of the periodic glitches and the seafloor noises with levels similar to the secondary microseism and higher than the NHNM (in the 3 × 10 -3 Hz -4 × 10 -2 Hz band) with a peak amplitude at 5 × 10 -3 Hz; Black: PSD of the Z-component after correction for the 1-hr glitches (black signal in supplementary Fig. S1a); The -130 dB peak at 5 × 10 -3 Hz remains, pointing to its oceanic origin; Its amplitude is weaker than the secondary microseism but still stronger than the primary microseism; Green: PSD of the Z-component after correction for the glitches and the tilt noise; Purple: PSD of the Z-component after correction for the glitches and seafloor noises;

The maximum amplitude of the residual signal is now well below the primary microseism and the NHNM. The dotted and dashed black curves are the PSDs of the horizontal components. Colored areas show the reduction of the noise level after correcting from glitches (blue), from the tilt (gray), and from the compliance (green). (b) Coherence between the vertical channel and the horizontal (blue and black curves) and pressure (red curve) channels. The gray area shows the frequency domain where the Z-component is coherent with the horizontal components due to tilt noise. The green area shows the domain where the Z-component is coherent with the pressure component due to compliance noise. At frequencies higher than 5 × 10 -2 Hz, the coherence is due to primary and secondary microseisms.

Accepted Article

This article is protected by copyright. All rights reserved. The Rayleigh waves have an average signal-to-noise ratio (SNR) greater than 3.5 (Supplementary Table S1). Although correlations for OBS-land station pairs are computed from only 8 months of data, Rayleigh waves have a SNR > 3 in the 5 -70 s period band, except in the 40 -70 s band where the SNR is slightly lower (SNR = 2.87, see supplementary Fig. S4b and Table S1). These noise correlations can therefore be used for Rayleigh-wave tomography. Supplementary Figures S4b-c demonstrate the effectiveness of corrections for glitches and seafloor noise in enhancing the SNR of correlations between OBS and land-station records.

First-order correlations between raw OBS records are displayed in the Fig. 3a. They can be compared to correlations for land-land station pairs and OBS-land station pairs shown in Supplementary Fig. S4. The SNR of correlation signals for OBS pairs is poorer than for other types of pairs in all period bands (average SNR < 2.6, Table 1). Therefore, Rayleigh waves are hardly detectable at periods shorter than 40 s in OBS-OBS correlations, and undetectable at longer periods (SNR < 1.6, Table 1). This may be explained by several factors such as local noises generated around the sensors by seafloor currents, or seismic noises generated between the stations that induce signals around time 0 s of the correlations, masking the Rayleigh waves. Therefore, these first-order OBS-OBS correlations of raw records cannot be used for Rayleigh-wave tomography.

Figure 3b and Table 1 document once again the effectiveness of the corrections for glitches and seafloor noises applied to OBS records, which improve the inter-OBS correlation signals. The correlations have a better SNR and are more symmetrical than those obtained from uncorrected signals, particularly at periods shorter than 40 s (columns a-b of Fig. 3).

However, correlation signals are still noisy, particularly at short lag times, and accurate measurements of Rayleigh-wave group velocities remain challenging. To further improve the quality of correlations for OBS pairs, we chose to virtually reconstruct the Rayleigh waves by computing iterative correlations.

Accepted Article

This article is protected by copyright. All rights reserved. Rayleigh waves emitted by each virtual source and recorded by two OBSs (i.e., by computing 272 a second-order correlation), it is possible, thanks to the stationary phase theorem, to isolate 273 the Rayleigh-wave propagating between the two OBSs, and therefore to measure its travel 274 time. We explain this with more detail in the following.

275

Let us consider any medium with a distribution of sources f . The wavefield recorded at a station A can be expressed using the Green's function G of the medium:

u( ⃗ r A , t) = Ω ∞ 0 G( ⃗ r s , ⃗ r A , t ′ ) f ( ⃗ r s , t -t ′ ) dt ′ d ⃗ r s (1)
It has been shown that the time-derivative of the first-order cross-correlation C 1 ⃗ r A , ⃗ r B (τ ) computed between wavefields recorded at two stations A and B is the Green's function

G ⃗ r A , ⃗ r B (t)
of the medium, assuming for instance a perfectly homogeneous distribution of white noise everywhere in the medium (e.g., [START_REF] Lobkis | On the emergence of the green's function in the correlations of a diffuse field[END_REF][START_REF] Snieder | Extracting the green's function from the correlation of coda waves: A derivation based on stationary phase[END_REF][START_REF] Wapenaar | Retrieving the elastodynamic green's function of an arbitrary inhomogeneous medium by cross correlation[END_REF][START_REF] Roux | P-waves from cross-correlation of seismic noise[END_REF][START_REF] Weaver | Information from seismic noise[END_REF][START_REF] De Verdière | Mathematical models for passive imaging i: general background[END_REF]:

d dτ C 1 ⃗ r A , ⃗ r B (τ ) = G ⃗ r A , ⃗ r B (τ ) (2)
In the case where all white noise sources are spatially uncorrelated, the correlation of wavefields recorded at A and B can be rewritten as the integral of correlations between G(⃗ r, ⃗ r A )

and G(⃗ r, ⃗ r B ):

C 1 ( ⃗ r A , ⃗ r B , τ ) = ∞ 0 u( ⃗ r A , t) u( ⃗ r B , t + τ ) dt = Ω ∞ 0 G( ⃗ r s , ⃗ r A , t ′ ) f ( ⃗ r s , t -t ′ ) dt ′ d ⃗ r s ⊗ Ω ∞ 0 G( ⃗ r s , ⃗ r B , t ′ ) f ( ⃗ r s , t -t ′ ) dt ′ d ⃗ r s = Ω G( ⃗ r s , ⃗ r A , t) ⊗ G( ⃗ r s , ⃗ r B , t) d ⃗ r s (3)
where ⊗ denotes the cross-correlation operation. Since the time-derivative of correlations is similar to the Green's function of the medium, it follows immediately by substituting dC 1 /dt

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Journal of Geophysical Research: Solid Earth to G in Eq. 3 that the first-order correlation is equivalent to a second-order correlation that we will note C 2 :

C 1 ( ⃗ r A , ⃗ r B , τ ) = Ω d dt C 1 ( ⃗ r s , ⃗ r A , t) ⊗ d dt C 1 ( ⃗ r s , ⃗ r B , t) d ⃗ r s = C 2 ( ⃗ r A , ⃗ r B , τ ) (4)
Equation 4 indicates that it is possible to reconstruct the Green's function of the medium between A and B by re-correlating the noise correlations computed between each point of the medium and stations A and B. However, this demonstration assumes that the timederivative of correlations are the exact and complete Green's function of the medium. The assumption would be correct if stations that could be used as virtual sources would exist everywhere in the medium.

In practice, the seismic noise recorded at the Earth surface in the period band considered in this work (5 -150 s) is dominated by Rayleigh waves. As a consequence, noise correlations computed between onshore stations and OBSs do not provide the full Green's function of the medium including all propagating modes, but they do provide robust estimates of the travel time of Rayleigh waves.

Rather than attempting to reconstruct the full Green's function, we will measure the travel time of the Rayleigh waves between OBSs A and B by using a simplified approach inspired by Eq. 4. It consists in re-correlating only the Rayleigh-wave parts of the correlations computed between onshore stations and OBSs A and B, that we consider as the Rayleigh-wave part of the Green's function:

d dτ C 2,Ray ( ⃗ r A , ⃗ r B , τ ) = d dτ Ω d dt C 1,Ray ( ⃗ r s , ⃗ r A , t) ⊗ d dt C 1,Ray ( ⃗ r s , ⃗ r B , t) d ⃗ r s ≈ G Ray ( ⃗ r A , ⃗ r B , τ ) (5)
where superscript "Ray" indicates that we only correlate the fundamental mode of Rayleigh waves of first-order correlations. In that way, we only retrieve the fundamental mode of the Rayleigh-wave part of the Green's function.

In practice, the distribution of virtual sources is never homogeneous. Instead, we use land stations deployed all over Western Europe, while OBSs are located in the Ligurian Sea (Fig. 1). Therefore, the condition of the stationary phase theorem are not completely met. To circumvent this difficulty, we select virtual sources that are expected to contribute constructively to the correlations (Fig. 4). To that end, we design a virtual-source azimuthalfilter that only retains sources located in the end-fire lobe of the OBS couple, that is in azimuths at ± 20 • with respect to the azimuth of the OBS pair. Moreover, we enhance the

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Journal of Geophysical Research: Solid Earth virtual source coverage by using separately the causal and anticausal parts of the first-order correlations to compute the second-order correlation.

Because onshore stations are well distributed and OBS-onshore correlations exhibit clear Rayleigh waves, the use of virtual sources in iterative correlations for OBS-OBS paths leads to a higher quality of Rayleigh waveforms (columns b-c of Fig. 3 and Table 1). This is achieved through: (1) separately recovering the causal and anticausal parts of the Green's function by using separately the causal and anticausal parts of the first-order correlations, thus avoiding interferences at long periods, and (2) controlling the distribution of virtual sources, thus guaranteeing a higher quality of the C 2 s by contrast to OBS-OBS C 1 s that exhibit low signal-to-noise ratio probably due to local noise sources. These local noise sources do not contribute significantly to OBS-onshore stations paths, and therefore neither to second-order correlations. Supplementary Fig. S5 shows that phases of the Rayleigh waveforms reconstructed from the C 1 and C 2 processes match. In the 5 -10 s band, iterative correlations do not systematically improve the signal quality as compared to firstorder correlations. Therefore, we select for each path the correlation of highest quality after checking that the C 1 and C 2 are coherent.

The strengths of the iterative correlations make it possible to substantially improve the path coverage in the Ligurian-Provence domain (Table 1). Our results demonstrate the efficiency of this method in providing robust group-velocity measurements. Further illustration and validation will be the subject of a future paper.

Group-velocity measurements

Once first-order and iterative correlations have been computed for onshore and offshore stations respectively, we derive group-velocity dispersion curves of the causal and acausal parts of the correlations by using multiple filter analysis (MFA, [START_REF] Dziewonski | A technique for the analysis of transient seismic signals[END_REF][START_REF] Herrmann | Some aspects of band-pass filtering of surface waves[END_REF]. As in [START_REF] Nouibat | Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps[END_REF], we adapt the width of the gaussian filter to the inter-station distance to accommodate the trade-off in resolution between the time and frequency domains [START_REF] Nature ; Levshin | Seismic Surface Waves in a Laterally Inhomogeneous Earth[END_REF]. We correct our group-velocity measurements from the biases that occur when the MFA method is applied on signals having a non-flat spectrum [START_REF] Shapiro | A systematic error in estimating surface-wave group-velocity dispersion curves and a procedure for its correction[END_REF]. This is especially important when measuring Rayleighwave velocities using noise correlations around the first and second microseismic peak, i.e., around 7 s and 14 s.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved. In order to build the group-velocity maps in the Ligurian-Provence domain, we maximize the path coverage over the Ligurian Sea by using simultaneously OBS-OBS, land-land, and land-OBS station pairs. A careful selection of group-velocity measurements is achieved to keep the most reliable ones and discard those that are biased by an unfavorable distribution of noise sources, or by interferences of causal and acausal Rayleigh waves for instance.

For first-order correlations (C 1 s) computed between land-land and land-OBS stations, at each period, we keep measurements if: (1) the SNR defined as the ratio of the Raleighwave peak amplitude and the standard deviation of the following signal, is greater than 3 on the positive and negative correlation times, (2) group velocities measured in positive and negative correlation times differ by less than 0.2 km/s, and (3) the inter-station distance is greater than 2 wavelengths. For iterative correlations (C 2 s) computed between OBS stations, we do not use the SNR criteria, since iterative correlations only exhibit the fundamental mode of the Rayleigh waves, owing to their construction. Nevertheless, for each C 2 satisfying the distance and symmetry criteria, we only keep the group-velocity measurement of the side of the correlation (positive or negative time) where the amplitude of the Rayleigh-wave is maximum (i.e., where we have more virtual sources contributing).

Table 1 shows that for OBS-OBS pairs, the selection procedure would reject more than 97% of Rayleigh waves velocity measurements performed on first-order correlations computed using OBS data that were not corrected from the compliance and tilt noises. This column (Supplementary Figure S7). The 10 s map still shows velocities lower than 2.5 km/s in the central and southwestern Ligurian-Provence basin. The 12 s and 15 s maps (sensitive to ∼ 10 -15 km depths) highlight velocities larger than 3.5 km/s associated with the thin crust of the Ligurian basin. However, the northeastern basin has lower velocities indicative

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 Journal of Geophysical Research: Solid Earth we will only summarize it here. We will rather focus on the specificity of the inversion for the offshore region, which is the consideration of the water layer. A result of the inversion for V s at an offshore location in the Ligurian Sea is shown in Supplementary Fig. S8.

Firstly, a 3-D probabilistic solution is computed that gives at each location the probability distribution of V s and the probability of having an interface as a function of depth. This is achieved using an exhaustive grid search on a set of ∼ 130 million synthetic four-layer models, which include a sedimentary layer, the upper crust, the lower crust and a half-space representing the upper mantle. The strength of this first-step Bayesian framework lies in constraining the structural complexity of the crust (i.e., of the short-period part of the dispersion curve) by means of an ensemble of models, fitting the dispersion curve to the degree required by its uncertainties. However, due to the four-layer model assumption, this procedure is not sufficient to fully describe the complexity of the model structure, particularly in the mantle part, that is in the long-period part of the dispersion curve. Hence, we use an iterative linear least-square inversion [START_REF] Herrmann | Computer Programs in Seismology: An Evolving Tool for Instruction and Research[END_REF] as a complement to update the mantle part of the model and further refine the fit to the local dispersion curve for the crustal part (Supplementary Fig. S8).

The initial model for the linear inversion (second step) at a given location is the average of all selected probabilistic solutions of the first inversion step, weighted by their likelihood values. We discretize the crustal and mantle parts at intervals of 1 and 5 km, respectively, and assume a gradual increase of V s below Moho according to the global model PREM [START_REF] Dziewonski | Preliminary reference earth model[END_REF]. For offshore locations, we incorporate on top of this initial model an additional layer of thickness equal to water depth at the given location. The parameters of this water layer are kept fixed during the linear least-square inversion (thickness,

V s = 0 km/s, V p = 1.5 km/s, ρ = 1 × 10 3 kg/m 3 ). Since we invert for short periods as well, an appropriate parameterisation of the water column is crucial. Indeed, we highlight its influence by computing group-velocity dispersion curves for synthetic 3-layer crustal models representative of the oceanic crust of the Ligurian-Provence basin, and different water levels from 0 to 3.1 km (Supplementary Fig. S7). This shows that the effect of the water layer on the dispersion curve is substantial at periods shorter than 15 -20 s, where the water depth changes impact: (1) the absolute group velocities, and (2) the shape of the dispersion curve, particularly in the vicinity of the Airy phase.

Accepted Article

This article is protected by copyright. All rights reserved. We think that such strong differences between the two V s models are mostly due to differences in coverage and quality of Rayleigh-wave dispersion data. Indeed, we greatly enhanced the path coverage in the Ligurian basin by using iterative correlations for OBS-OBS paths. In addition, the use of all broadband stations in Western Europe provides long ray paths across the Ligurian sea, allowing us to improve path coverage to the west and southeast and make use of Rayleigh-wave dispersion measurements up to 150 s. On the other hand, Wolf et al. ( 2021) used only 23 onshore stations nearby the basin, thus a limited aperture. Therefore, they used earthquake records for periods greater than 20 s.

To a lesser extent, the differences between the two models may also be explained by the different strategies used to invert Rayleigh-wave dispersion measurements. Wolf et al. (2021) computed their 2-D dispersion maps using a linear inversion that depends on an explicit regularization, while our transdimensional approach does not. Moreover, our 1-D Bayesian inversion for V s take uncertainties on dispersion measurements into account, which is key for controlling model complexity.

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Journal of Geophysical Research: Solid Earth

As we will see in the following section, our V s model is more coherent with current knowledge on the crustal structure of the Ligurian-Provence basin than the one by Wolf et al. (2021). In particular, the thickness of the sedimentary cover and Moho depth and geometry estimated from our V s model are coherent with the V p model by [START_REF] Dannowski | Seismic evidence for failed rifting in the ligurian basin, western alpine domain[END_REF] along the basin axis and with the stratigraphic log after [START_REF] Leprêtre | Multiphased tectonic evolution of the central algerian margin from combined wide-angle and reflection seismic data off tipaza, algeria[END_REF].

Discussion

We now focus on the oceanic domain of the Ligurian-Provence basin (Fig. 8a), in particular on the comparison of our S-wave velocity model with a recent P -wave velocity model derived by [START_REF] Dannowski | Seismic evidence for failed rifting in the ligurian basin, western alpine domain[END_REF] from a controlled-source seismic profile recorded along the basin axis (thick black line in Fig. 8a). The availability of this high-resolution

V p section provides a unique opportunity to assess the accuracy and validate our V s model against an independent dataset. Moreover, the existence of V p and V s models along the same profile may provide clues on the petrological structure of the crust in the oceanic domain of the basin, which is still debated. Indeed, it has been proposed that the oceanic domain is made of an oceanic crust with a thin basaltic layer (e.g., [START_REF] Mascle | A revised seismic stratigraphy of the tyrrhenian sea: implications for the basin evolution[END_REF][START_REF] Bonatti | Peridotites drilled from the tyrrhenian sea, odp leg 107[END_REF], or an exhumed and serpentinized mantle devoid of any volcanic upper layer (e.g., [START_REF] Boillot | Undercrusting by serpentinite beneath rifted margins[END_REF][START_REF] Beslier | Ocean-continent boundary in the iberia abyssal plain from multichannel seismic data[END_REF]Jolivet et al., 2020), or even an hyper-extended continental crust (e.g., [START_REF] Mckenzie | Some remarks on the development of sedimentary basins[END_REF][START_REF] Pascal | The ocean-continent boundary in the gulf of lion from analysis of expanding spread profiles and gravity modelling[END_REF][START_REF] Dannowski | Seismic evidence for failed rifting in the ligurian basin, western alpine domain[END_REF].

Geological setting of the Ligurian-Provence basin

The basin opening initiated at 30 Ma by a rifting phase between Europe and the Corsica-Sardinia block, as a result of back-arc extension above the Adria oceanic microplate, initially subducting north-westward (e.g., [START_REF] Faccenna | Styles of back-arc extension in the central mediterranean[END_REF]. The progressive south-eastward roll-back and retreat of the Adria slab below the Corsica-Sardinia domain led to stretching of the continental crust followed by continental break-up during the early Miocene, and to the genesis of an oceanic crust between 20 and 15 Ma [START_REF] Séranne | The gulf of lion continental margin (nw mediterranean) revisited by ibs: an overview[END_REF].

As a result, the Ligurian-Provence basin includes two thinned conjugate continental passive margins separated by an oceanic domain (Fig. 8a). According to [START_REF] Rollet | Back arc extension, tectonic inheritance, and volcanism in the ligurian sea, western mediterranean[END_REF], the entire region is characterized by magnetic anomalies, and by the presence of magmatic bodies identified from acoustic facies in seismic reflection profiles. The area between the margins and the oceanic domain is described as a transitional domain, likely made up of a very thin continental crust overlying a thick rift-related corner of magmatic underplating (e.g.,

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 Journal of Geophysical Research: Solid Earth gabbroic intrusion (Fig. 9c). Owing to the presence of the gabbro intrusion, a single S-wave velocity contour cannot be used as proxy for the petrological Moho in the Ligurian-Provence basin, unlike in continental areas [START_REF] Nouibat | Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps[END_REF].

In the southwestern part of the profile, the depths of the V s (4.1 km/s) and V p (7.2 km/s) Moho proxies differ by less than 1 km. Such a small discrepancy is remarkable, given that the two models are totally independent. The depth profiles of the two shallower layer boundaries are also remarkably similar to those of the V p contours that define lithological layering in the western Mediterranean basins (Fig. 8b). Such similarity to the P -wave velocity model of [START_REF] Dannowski | Seismic evidence for failed rifting in the ligurian basin, western alpine domain[END_REF] validates the offshore part of our shear-wave velocity model, as similarity to the receiver function section of the Cifalps profile validated its onshore part [START_REF] Nouibat | Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps[END_REF].

In the few locations where information on P -wave velocity is available, the uppermost mantle has the seismic signature of a dry peridotite, with V p > 7.2 km/s and V s > 4.0 km/s [START_REF] Grevemeyer | Episodic magmatism and serpentinized mantle exhumation at an ultraslow-spreading centre[END_REF]. We find no evidence of serpentinized mantle, which would show much lower P -and S-wave velocities. [START_REF] Dannowski | Seismic evidence for failed rifting in the ligurian basin, western alpine domain[END_REF] and Wolf et al. ( 2021) also concluded on low mantle serpentinization in this part of the basin axis.

Conclusion

Using data of 23 ocean-bottom seismometers of the AlpArray network with those of 890 temporary and permanent onshore stations, we have derived a 3-D high-resolution shearwave velocity model encompassing the Ligurian-Provence basin and its conjugate margins.

The OBS continuous records could be fully exploited after a careful, specific pre-processing scheme including removal of instrument noises (glitches) and reduction of seabed-induced compliance and tilt noises. We enhanced the quality of correlations between OBSs and maximized the path coverage in the Ligurian-Provence basin by involving correlations with onshore stations to virtually reconstruct Rayleigh waves propagating between OBSs. As in [START_REF] Nouibat | Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps[END_REF], we computed 2-D group-velocity maps and their uncertainties using a data-driven transdimensional inversion of Rayleigh-wave group-velocity measurements. The dispersion data and their uncertainties have then been used jointly in a Bayesian probabilistic approach to derive a 3-D probabilistic shear-wave velocity model. The output average model was further refined using a linear inversion that accounts for the presence of the water column.

Accepted Article

This article is protected by copyright. All rights reserved. The comparison with the high-resolution P -wave velocity section derived by [START_REF] Dannowski | Seismic evidence for failed rifting in the ligurian basin, western alpine domain[END_REF] from traveltime inversion of controlled-source seismic data along the basin axis validates our 3-D ANT model. Layer boundaries revealed by high probabilities of presence of an interface and V s contours are remarkably consistent with V p contours. The joint interpretation of the V p and V s models highlights a relatively thin anomalous oceanic crust of low P -wave velocities but rather high S-wave velocities. In the NE part of the profile, the lower part of the crust exhibits a gabbroic intrusive body. The underlying mantle is anhydrous and shows no evidence of serpentinisation. These results show the potential of a joint interpretation of V p and V s models since they provide reliable answers to a number of debated questions on the petrological nature of the crust and uppermost mantle of the Ligurian-Provence basin, at least along the LOBSTER-P02 seismic profile.

They also warrant the same type of study on the SEFASILS controlled-source seismic profile that crosses the northern margin of the basin [START_REF] Dessa | Seismic exploration of the deep structure and seismogenic faults in the ligurian sea by joint multi channel and ocean bottom seismic acquisitions: Preliminary results of the sefasils cruise[END_REF].

The use of OBS recordings in ambient-noise tomography is more challenging than with onshore stations due to shorter recording times, a higher potential of technical problems, sea-floor noises and generally a poorer signal-to-noise ratio in the frequency bands of microseismic noise that are key for ambient-noise tomography. Parts of these problems have been solved here by a specific pre-processing of OBS records that reduces instrument and sea-floor generated noises such as tilt and compliance. When available, land stations can be used in combination with OBSs to provide higher quality surface-wave signals between offshore and onshore stations than for OBS pairs, therefore improving interstation path coverage of ANT in particular at the ocean-continent transition. We went a step further by showing how to take full benefit of onshore stations to enhance the quality of correlations for OBS pairs. Our Rayleigh-wave reconstruction scheme for OBS pairs based on second-order correlations between OBSs and land stations has proven to be effective in improving the coverage of the offshore domain. Finally, the transdimensional inversion of the enhanced set of Rayleighwave group-velocity observations for group-velocity maps and their uncertainties and the following hybrid inversion for V s that accounts for the water layer have lead to a high-quality 3-D V s model of the study region. We have therefore set up a complete, efficient and reliable ambient-noise imaging methodology of oceanic domains and their margins using OBSs and land stations that opens new perspectives for the processing of similar datasets.

  manuscript submitted to Journal of Geophysical Research: Solid Earth

  Figure 1. Tectonic map of the western Mediterranean region (modified from Faccenna et al., 2014; Jolivet et al., 2020) with locations of seismic stations used in this work (red circles: onshore sensors; green triangles: AASN OBS). The oceanic domains of the Ligurian-Provence, Algerian and Tyrrhenian back-arc basins are filled with gray color. The gray frame shows the location of the map in Figure 8a.

  manuscript submitted to Journal of Geophysical Research: Solid Earth 2 Data processing 2.1 Description of the AASN sea-bottom instruments

  and Earthquake Research) belong to the GEOMAR and DEPAS pools (German instrument pool for amphibious seismology). Designed by K.U.M. Enviromental-and Marine Technology GmbH, these instruments were equipped with HTI-01-PCA hydrophones from High Tech Inc., Trillium compact velocimeters from Nanometrics and K.U.M., and with CMG-40T velocimeters from Güralp for four of them. Seven broadband OBS (BBOBS) belong to the French OBS pool of INSU-IPGP. Designed by the Scripps Institution of Oceanography, these instruments were equipped with deep-sea differential pressure gauges and Nanometrics Trillium-240 very broadband velocimeters.

Figure 2 .

 2 Figure 2. Effect of the different corrections for glitch and seafloor noises on the PSD of a 1day vertical-component record of OBS A416A. Raw signals are band-pass filtered between 2.5 and 250 s, corrected from the instrumental response and decimated to 1 Hz sampling frequency. (a) PSDs before and after corrections; Orange dashed lines indicate the Peterson (1993) New Low and High Noise Model (NLNM & NHNM); Blue: PSD of the raw vertical-component record (blue signal in supplementary Fig. S1a); At frequencies > 5 × 10 -2 Hz, the PSD is dominated by the primary

  manuscript submitted to Journal of Geophysical Research: Solid Earth single reference correlation per station pair. We use up to 4 years of continuous vertical records for on-land stations pairs, and up to 8 months data to compute correlations for OBS-OBS and OBS-land-station pairs. As shown by[START_REF] Nouibat | Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps[END_REF] and supplementary Fig.S4a, Rayleigh waves are clearly visible in the correlations for on-land station pairs in a wide period band (5 -150 s).

  manuscript submitted to Journal of Geophysical Research: Solid Earth3.2 Iterative correlations for OBS-OBS paths 264We have seen in the previous section that correlations for OBS pairs can hardly be 265 used for Rayleigh-wave tomography. Since correlations computed between onshore stations 266 and OBSs exhibit clear Rayleigh waves, we will use onshore stations as virtual sources 267 in order to measure the travel time of Rayleigh waves between OBSs. Indeed, first-order 268 correlations computed between onshore stations and OBSs contain the Rayleigh-wave part 269 of the Green's function. It is thus possible to use them to mimic the case where Rayleigh 270 waves are emitted on the continent and recorded by OBSs. By computing correlations of 271

Figure 3 .

 3 Figure 3. Time-distance plots of correlation signals for OBS-OBS pairs, obtained in different period bands at different steps of the processing. (a) First-order cross-correlations (C 1 s) of raw signals. (b) C 1 s of pre-processed signals (glitch removal & seafloor-noise reduction). (c) Secondorder correlations (C 2 s). The C 1 s in (a) are generally of poor quality and poor symmetry, with strong signals at short lag times due to interferences even at short periods. These interferences become stronger with increasing period, and they progressively overshadow the Rayleigh wavetrains. In (b), the C 1 signals are strongly improved due to corrections for glitches and seafloor noises, in particular in the three short-and medium-period bands. The 70 -150 s band is still affected by interferences that hide the Rayleigh wave-trains. In (c), the quality and the timesymmetry of correlation signals are significantly improved in all period bands by the calculation of C 2 s. In particular, the Rayleigh wave-train emerges from the noise in the long-period band.

Figure 4 .

 4 Figure 4. Details of the computation of second-order correlations for 2 examples of OBS-OBS pairs, using station CFF (FR network) as virtual source. (1-2a): Station location maps. In example (1), the virtual source is roughly aligned with the OBS pair, while in example (2), the azimuth of the virtual source is almost perpendicular to the pair. (1-2b, 1-2c): First-order correlation signals between CFF and each OBS. The green areas show the Rayleigh-wave search windows (wave propagation of 1 to 5 km/s over the inter-station distances). The blue areas show the Rayleigh-wave detection windows. (1-2d): First-order correlation of the OBS signals (C 1 , in black), and correlation of the positive-time (causal) parts of the C 1 s between CFF and the OBS couple (P-P, in blue).(1-2e): First-order correlation of the OBS signals (C 1 , black), and correlation of the negative-time (acausal) parts of the C 1 s between CFF and the OBS couple (N-N, red). (1-2f): Comparison of the causal (P-P, blue) and acausal (N-N, red) C 2 s (that will be summed to obtain the final C 2 ), with the first-order correlation of the OBS signals (C 1 , black). The C 2 and C 1 signals have similar phases in (1) where the virtual source is aligned with the OBS pair. They have incoherent phases in (2), since CFF is not located in the end-fire lobe of this OBS pair (see text).

Figure 5 .

 5 Figure 5. Group-velocity maps (average solutions) at 6 to 35 s periods, obtained with the Hierarchical Bayes reversible-jump algorithm. Only areas with uncertainty lower than 0.5 km/s are shown. C-SLPB: central and southwestern parts of the basin, GL: Gulf of Lion, NLPB: northeastern Ligurian-Provence basin, SFB: southeast-France basin.

  for shear-wave velocity 383 The group-velocity maps and their uncertainties are used to derive a 3-D V s model. For 384 this, we perform a two-step data-driven inversion to tackle the non-unicity of the inverse 385 problem. The main part of this process is described in detail in Nouibat et al. (2022), so 386 -19-

Figure 6 .

 6 Figure 6. Depth slices in the final Vs model at 3 to 30 km depths. Only regions with 1σ error < 8 % are shown.

Figure 6

 6 Figure6shows depth slices in the 3-D shear-wave velocity model at 3 to 30 km depth.

418

  The sediment layer is clearly visible at 3 -7 km depth in the central and southwestern 419 Ligurian-Provence basin, with low velocities of 2 -3.2 km/s. The transition from crustal 420 (3.5 -4.1 km/s) to mantle velocities (4.1 -4.5 km/s) is located between 10 km and 12 -421 15 km in the parts of the basin with thinnest crust. Between 20 and 30 km, areas of strong 422 velocities corresponding to the mantle extend westward from the Ligurian-Provence basin 423 to the Gulf of Lion across a nearly N-S transition of slower velocities at ∼ 5.5 • E.

4244. 3 Figure 7

 37 Figure 7 shows a comparison of our V s model and the model by Wolf et al. (2021),

  manuscript submitted to Journal of Geophysical Research: Solid Earth the Ligurian basin. Our V s model covers a wider area as part of the large-scale model by[START_REF] Nouibat | Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps[END_REF] that uses all available broadband stations in Western Europe.At 3-km depth, our model (Fig.7-1a) highlights: (1) the sedimentary cover in the central and southwestern basin (1.8 -2.2 km/s), where Wolf et al's model (Fig. 7-1b) exhibits patches of very low velocities (V s ≤ 1.4 km/s), and (2) gradual increasing of velocity from the central basin towards the conjugate margins, which is not clearly visible in Wolf et al's model where velocities fluctuate from very low (V s ≤ 1.8 km/s) to high values (2.4 -2.5 km/s). At 5-km depth, we still observe typical sediment velocities (Fig. 7-2a), while Wolf et al's model (Fig. 7-2b) exhibits higher velocities (V s > 3.5 km/s). At 13-km depth, our model (Fig. 7-3a) shows almost homogeneous mantle velocities in the basin (V s ≥ 4.1 km/s) while Wolf et al's model (Fig. 7-3b) exhibits numerous small-size velocity anomalies, fluctuating between mantle-like (V s ≥ 4.1 km/s) and crust-like velocities (V s ≤ 3.7 km/s). These heterogeneities suggest a much more irregular Moho surface than in our model. At 20-km depth, our model (Fig. 7-4a) documents crustal thinning towards the basin axis with V s increasing across the conjugate margins. Again, velocities in the central basin are more heterogeneous in Wolf et al's model (Fig. 7-4b), with localized small-size anomalies. Although the transition domain has an irregular shape in Wolf et al's model, it remains to first order similar to that shown in our model.

Figure 7 .

 7 Figure 7. Comparison between our Vs model (left panel) and the Vs model by Wolf et al. (2021) (right panel) at four depths.

  manuscript submitted to Journal of Geophysical Research: Solid Earth

Table 1 .

 1 Signal-to-noise ratio and percentage of selected paths for group-velocity tomography

	manuscript submitted to Journal of Geophysical Research: Solid Earth
		SN R -% of retained paths
	Period band 1st order CC Cleaned 1st order CC 2nd order CC
	10s -20s	2.52 -2.45%	3.29 -8.11%	5.77 -39.6%
	20s -40s	2.32 -1.52%	3.61 -11.5%	5.62 -23.4%
	40s -70s	1.57 -0.03%	3.12 -5.22%	4.52 -14.3%
	70s -150s	1.13 -0.01%	1.42 -2.61%	4.37 -7.31%
	in different period bands using: (1) first-order cross-correlations of raw OBS vertical-component
	records, (2) first-order cross-correlations of pre-processed signals (corrected for glitches, and seafloor
	noises), and (3) second-order cross-correlations of pre-processed signals.

This article is protected by copyright. All rights reserved.

Acknowledgments

This study was partly funded by the SEISCOPE consortium (https://seiscope2.osug.fr), sponsored by AKERBP, CGG, CHEVRON, EQUINOR, EXXON-MOBIL, JGI, SHELL, SINOPEC, SISPROBE and TOTAL, and by the RGF program (Référentiel Géologique de la France, (http://rgf.brgm.fr/). It is part of the AlpArray-FR project funded by Agence Nationale de la Recherche (contract ANR-15-CE31-0015) and by Labex OSUG@2020 (Investissements d'Avenir, ANR-10-LABX-56). The deployment of OBSs for 8 months in the Liguro-Provençal basin was funded by project AlpArray-FR for the French component, and by the LOBSTER project that is part of the German Priority Programme SPP2017 4D-MB for the German component. This study was granted access to the HPC resources of the Dahu platform of the CIMENT infrastructure (https://gricad.univ-grenoble-alpes.fr/), which is supported by the Auvergne-Rhône-Alpes region (grant CPER07-13CIRA), the Labex OSUG@2020 (reference ANR10 LABX56) and the Equip@Meso project (reference

Data availability statement

Our tomographic 3-D model will be available in the Résif products repository (Réseau sismologique et géodésique Français; https://www.resif.fr/en/data-and-products/productsrepository/). The inversion code is available on request to the first author.

Waveform data used in this paper are available through the EIDA (European Integrated Data Archive) service of ORFEUS (http://www.orfeus-eu.org/eida/) and belong to the permanent networks with codes AC, BE (Royal Observatory Of Belgium, 1985), CA (Institut Cartogràfic I Geològic De Catalunya -Institut D'Estudis Catalans, 1984), CH (Swiss Seismological Service (SED) At ETH Zurich, 1983), CR (University Of Zagreb, 2001), CZ (Institute Of Geophysics Of The Academy Of Sciences Of The Czech Republic, 1973), ES (Instituto Geografico Nacional, Spain, 1999), FR (RESIF, 1995), G (Institut de Physique du Globe de Paris (IPGP) & Ecole et Observatoire des Sciences de la Terre de Strasbourg (EOST), 1982), GB (GEOFON Data Centre, 1993), GR (Federal Institute for Geosciences & Natural Resources (BGR), 1976), IV (INGV Seismological Data Centre, 2006), NL (KNMI, 1993), OE (ZAMG-Zentralanstalt Für Meterologie Und Geodynamik, 1987), SL (Slovenian Environment Agency, 2001), and UP (SNSN, 1904). We also used data of the temporary Al-pArray network (network code Z3, AlpArray Seismic Network, 2015), Cifalps-2 experiments (network code XT, Zhao et al., 2018) and EASI experiments (network code XT, AlpArray Seismic Network, 2014).

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Journal of Geophysical Research: Solid Earth illustrates that these signals are obviously not usable for ambient-noise tomography. By contrast, we selected between 2.6% to 11.5% of the measurements performed on first-order correlations done using OBS data corrected from sea floor noises. This highlights the importance of the pre-processing scheme described in section 2. Moreover, we kept between 7.3% and 39.6% of group-velocity measurements performed on iterative correlations (C 2 s) depending on the period-band considered. Second-order correlations provide a substantial gain over first-order correlations, leading to a significant improvement of the raypath coverage in the Ligurian-Provence basin with respect to the ANT of Wolf et al. (2021) based on first-order correlations.

4 3-D shear-wave velocity model

Inversion for 2-D group-velocity maps

We compute 2-D group-velocity maps and associated uncertainties using a "datadriven" transdimensional approach at discrete periods from 4 s to 150 s. At each period, probabilistic group-velocity maps are derived by exploring millions of 2-D models using the reversible-jump Markov-chain Monte-Carlo method (rj-McMC, [START_REF] Bodin | Transdimensional tomography with unknown data noise[END_REF]. The method used for the inversion and the spatial resolution of the resulting group-velocity maps are discussed in detail in [START_REF] Nouibat | Lithospheric transdimensional ambient-noise tomography of W-Europe: implications for crustal-scale geometry of the W-Alps[END_REF]. Uncertainty and path density maps at different periods are presented in Supplementary Figure S6.

Resulting group-velocity maps of the Ligurian-Provence basin and its margins are shown in Fig. 5 for periods from 6 s to 35 s. The 6 s and 8 s maps (sensitive to ∼ 4 -8 km depth) highlight low-velocity anomalies (U < 1.6 km/s) in the central and southwestern parts of the basin (C-SLPB in Fig. 5), that are probably associated with thick sediment sequences.

These velocities are lower than those of the southeast-France basin and the Gulf of Lion (respectively SFB and GL in Fig. 5). The northeastern Ligurian-Provence basin (NLPB), its northern, Provence coast margin and its southern margin in Corsica have larger velocities (U ≥ 2.4 km/s). From simulations of Rayleigh-wave dispersion in synthetic 1-D models, we show that such periods are also highly sensitive to the presence and thickness of the water

Accepted Article

This article is protected by copyright. All rights reserved.

manuscript submitted to Journal of Geophysical Research: Solid Earth [START_REF] Séranne | The gulf of lion continental margin (nw mediterranean) revisited by ibs: an overview[END_REF]. This limit is marked by an abrupt change in the amplitude of magnetic anomalies with a transition from mostly positive values in the deep basin (i.e., oceanic domain) to negative values at the continent-ocean transition, and by a change in acoustic facies on seismic reflection profiles (e.g., [START_REF] Réhault | The western mediterranean basin geological evolution[END_REF][START_REF] Déverchère | Malis cruise, rv le nadir[END_REF][START_REF] Rollet | Back arc extension, tectonic inheritance, and volcanism in the ligurian sea, western mediterranean[END_REF]. While the magmatism occurring in the margins has been associated to back-arc magmatic activity strongly influenced by subduction (e.g., [START_REF] Coulon | Le volcanisme calco-alcalin cénozoique de sardaigne (italie): Pétrologie, géochimie et genese des laves andésitiques et des ignimbrites: Signification géodynamique, these de doctorat[END_REF][START_REF] Bellon | Chronologie radiométrique (k-ar) des manifestations magmatiques autour de la méditerranée occidentale entre 33 et 1 ma[END_REF][START_REF] Réhault | Offshore oligo-miocene volcanic fields within the corsica-liguria basin: Magmatic diversity and slab evolution in the western mediterranean sea[END_REF], the nature of magmatism observed in the oceanic domain remains unknown. It could not be investigated by direct geochemical analysis due to the presence of a sedimentary cover several kilometres thick.

Figure 8b shows a stratigraphic log representative of Western Mediterranean oceanic basins that includes P -wave velocity estimates. It is derived from the results of joint seismic wide-angle and reflection profiling in the Algerian and Western Sardinia basins [START_REF] Klingelhoefer | Preliminary results from the sardinia deep seismic cruise on the western sardinia and gulf of lions conjugate margin pair[END_REF][START_REF] Gailler | Crustal structure of a young margin pair: New results across the liguro-provencal basin from wide-angle seismic tomography[END_REF][START_REF] Leprêtre | Multiphased tectonic evolution of the central algerian margin from combined wide-angle and reflection seismic data off tipaza, algeria[END_REF]. The sedimentary layer of 5-km average thickness and 1.9 -5 km/s P -wave velocities is made up of Plio-quaternary sediments, Messinian and pre-salt units (Fig. 8b). The transition from Plio-quaternary to Messinian units occurs at V p ≈ 2.5 km/s. The Messinian sequence exhibits strong thickness variations ascribed to salt diapirism. It is separated from the pre-salt unit by the 4.2 km/s velocity boundary. The deepest sediments overlay an oceanic basement that starts at V p > 5 km/s. The oceanic crust is relatively thin with an average thickness of 4.5 km. Its P -wave velocities range from 5 to 7.2 -7.3 km/s at ∼ 12 km depth, which corresponds to the Moho depth (V p > 7.3 km/s in the upper mantle). Controlled-source seismic data are useful to constrain the layer thicknesses and the depths of major interfaces (intra-sedimentary, sediment-crust and Moho), but their interpretation in terms of petrology only rely on Pwave velocity estimates. The interpretation is ambiguous as two lithologies of different petrological natures may have similar V p (or V s ) signatures. However, P -and S-wave velocities can be used jointly to yield information on lithologies and their hydration degree (e.g., [START_REF] Grevemeyer | Episodic magmatism and serpentinized mantle exhumation at an ultraslow-spreading centre[END_REF][START_REF] Malusà | The deep structure of the alps based on the cifalps seismic experiment: A synthesis[END_REF]. For instance, the V p /V s ratio is commonly used to assess the degree of serpentinization in oceanic domains and margins (e.g., [START_REF] Bullock | From continental extension to seafloor spreading: crustal structure of the Goban Spur rifted margin, southwest of the UK[END_REF]Reynard, 2013;[START_REF] Grevemeyer | Episodic magmatism and serpentinized mantle exhumation at an ultraslow-spreading centre[END_REF]. We will take advantage of the availability of the V p cross-section by [START_REF] Dannowski | Seismic evidence for failed rifting in the ligurian basin, western alpine domain[END_REF] and our V s section along the same profile to further constrain the petrological nature of the crust in the central Ligurian-Provence basin.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved. Figure 9b shows two major layer boundaries at ∼ 5 and ∼ 12 km depth with rather high probabilities of presence, and a third one of weaker probability at ∼ 7.5 km depth. The good correspondence of the shallowest boundary with the velocity contour V s = 2.5 km/s suggests that it is probably an intra-sedimentary interface. This boundary also coincides with the 4.2 km/s P -wave velocity contour (Fig. 9c), which corresponds to the base of the Messinian salt unit according to Fig. 8b. This interpretation is consistent with S-and Pwave velocities, typical of salt (e.g., Yan et al., 2016). The intermediate interface at ∼ 7.5 km depth, which is slightly less pronounced that the Moho boundary (Fig. 9b) coincides with the velocity contours V s = 3.5 km/s (Fig. 9b) and V p = 5 km/s (Fig. 9c), which may support its interpretation as the sediment-crust boundary, in agreement with Fig. 8b. However, this interface does not correspond to a marked change in seismic velocity in the final model of Fig. 9a. Therefore, we cannot detect unambiguously the depth of the sediment-crust transition. We will use the proxy V s = 3.5 km/s (or V p = 5 km/s) for the sediment-crust boundary.

In the northeastern part of the transect (x > 70 km), the lower crust has higher, but still crustal P -wave velocities (6.0 -7.2 km/s, in green in Fig. 9c), hence a weaker P -wave velocity gradient at the Moho than in the southwestern part. By contrast, S-wave velocities are high and almost mantle-like (4.0 -4.4 km/s), with a low S-wave velocity gradient at Moho depth, in particular at the northeastern end of the profile (x > 90 km in Fig. 9a).

These P -and S-wave velocities are typical of gabbro [START_REF] Grevemeyer | Episodic magmatism and serpentinized mantle exhumation at an ultraslow-spreading centre[END_REF], which suggests a gabbro intrusive body within the oceanic crust. Our interpretation is at odd with [START_REF] Dannowski | Seismic evidence for failed rifting in the ligurian basin, western alpine domain[END_REF] who interpreted the high V p part of the deep crust as hyperextended continental crust based on gravity modeling. The observed high S-wave velocity rules out the continental crust hypothesis.

As outlined by [START_REF] Dannowski | Seismic evidence for failed rifting in the ligurian basin, western alpine domain[END_REF], the V p = 7.2 km/s is a good Moho proxy because it coincides with a very strong velocity gradient. In the southwestern part of the profile (x < 70 km), the V p = 7.2 km/s contour closely corresponds to the V s = 4.1 km/s contour while it corresponds to the V s = 4.4 km/s contour for x > 70 km, that is beneath the

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article

This article is protected by copyright. All rights reserved.


Vasiliki Stergiopoulou

José Henrique De

M Goulart

Sébastien Schaub

Luca Calatroni

Laure Blanc-Féraud

José Henrique De Morais Goulart

COL0RME: COvariance-based

Keywords: ℓ_0 super-Resolution Microscopy with intensity Estimation Super-Resolution, Sparse Non-Convex Optimization, Fluorescence microscopy, SOFI, SRRF

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Fluorescence microscopy is an imaging technique that allows the investigation of living cells and their organelles. However, due to light diffraction, structures smaller than half the emission wavelength (typically, around 250 nm) cannot be resolved. As many entities of biological interest have a size smaller than such a barrier, it is therefore fundamental to consider approaches which are able to provide a super-resolved version of the acquired data. Techniques such as Single Molecule Localization Microscopy (SMLM) [START_REF] Sage | Quantitative evaluation of software packages for singlemolecule localization microscopy[END_REF] and STimulated Emission Depletion (STED) [START_REF] Hell | Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy[END_REF] can achieve nanometric spatial resolution at the cost of low temporal resolution and need of special (often costly) equipment. On the other side, Structured Illumination Microscopy (SIM) [START_REF] Gustafsson | Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy[END_REF] has fast acqui-sition rates but relatively low spatial acquisition.

To overcome these limitations, different super-resolution approaches taking advantage of the independent stochastic temporal fluctuations of standard fluorescent dyes were considered during the last decade. These methods take as an input a temporal stack of images acquired by common microscopes and produce a super-resolved image on a finer grid. In Super-resolution Optical Fluctuation Imaging (SOFI) [START_REF] Dertinger | Fast, background-free, 3D super-resolution optical fluctuation imaging (SOFI)[END_REF], for instance, fluctuations are exploited by mapping the acquired image data into the covariance domain. While SOFI can drastically improve the temporal resolution of the acquired images, its spatial resolution is still limited in comparison to SMLM and STED. Better results can be achieved by Super-Resolution Radial Fluctuations (SRRF) [START_REF] Gustafsson | Fast live-cell conventional fluorophore nanoscopy with ImageJ through super-resolution radial fluctuations[END_REF] microscopy, where a degree of local symmetry is computed for each frame. Recently, a SPARsity-based super-resolution COrrelation Microscopy (SPARCOM) [START_REF] Solomon | SPARCOM: Sparsity based super-resolution correlation microscopy[END_REF] approach has been proposed to exploit the lack of correlation between different emitters as well as the sparse distribution of the fluorescent molecules via the use of a sparsity-promoting 1 -type regularization of the emitters' autocovariance.

In the same spirit as SPARCOM, we propose a method for COvariance-based 0 super-Resolution Microscopy with intensity Estimation (COL0RME) where signal sparsity is enforced in the covariance domain. Differently from SPAR-COM, we use in our work an 0 -type penalization instead of the 1 norm of the signal temporal variance. Furthermore, a joint estimation of the noise variance is performed for a more precise modeling. Our approach further incorporates a second step which allows for both signal intensity and constant background estimation, in order to deal with realistic experimental conditions. The retrieval of real image intensity information is particularly relevant for biological interpretation. To the best of our knowledge, COL0RME is the only superresolution method exploiting temporal fluctuations which is capable of retrieving this information. Simulated and experimental results show that COL0RME performs well in terms of molecule localization and can retrieve accurate intensity information.

COL0RME: FORMULATION

Inverse problem formulation

Let Y t ∈ R M ×M be the image frame acquired at time t ∈ {1, . . . , T }, with T > 1 and X t ∈ R L×L , with L = qM , the desired high-resolution image defined on a q-times finer grid, with q ∈ N. The discrete model describing the acquisition process at frame t is given by:

Y t = M q (H(X t )) + B + N t (1)
where M q : R L×L → R M ×M is a down-sampling operator averaging every q consecutive pixels in both dimensions, H : R L×L → R L×L is a convolution operator defined by the PSF of the optical imaging system, B models the spatially and temporally constant background (due to out-of-focus fluorescent molecules) and N t denotes electronic noise, which is modeled here as a matrix of independent and identically distributed (i.i.d.) Gaussian elements, with constant variance s ∈ R + . We further assume that the molecules are located at the center of each pixel and that there is no displacement of the specimen during the imaging period.

In vectorized form the model (1) reads:

y t = Ψx t + b + n t (2)
where now Ψ ∈ R M 2 ×L 2 is the matrix representing the composition M q •H and lower-case letters imply the column-wise vectorization of the corresponding matrices in [START_REF] Sage | Quantitative evaluation of software packages for singlemolecule localization microscopy[END_REF]. Given the matrix Ψ describing the acquisition process and the frame y t , the task thus consists in computing a super-resolved image x from the many y t acquisitions and, if possible, in estimating b and s. In order to exploit the statistical behaviour of the fluorescent emitters we reformulate the model in the covariance domain in the next section.

Sparse regularization in the covariance domain

We exploit temporal and spatial independence of the fluorescent emitters by computing spatial statistics. This idea was previously exploited in [START_REF] Dertinger | Fast, background-free, 3D super-resolution optical fluctuation imaging (SOFI)[END_REF] and was shown to significantly reduce the full width at half maximum (FWHM) of the PSF. In particular, the use of second-order statistics corresponds to a reduction factor √ 2. We consider the frames y t as T realizations of a random variable y whose covariance matrix is defined as:

R y = E y {(y -E y {y})(y -E y {y}) } (3)
where E y {•} denotes the expected value computed w.r.t. to the unknown law of y. We approximate R y by its sample average:

R y ≈ 1 T -1 T t=1 (y t -y)(y t -y)
where y = 1 T T t=1 y t is the empirical temporal mean. From ( 2) and (3) we thus have:

R y = ΨR x Ψ + R n (4)
where

R x ∈ R L 2 ×L 2 and R n ∈ R M 2 ×M 2
are the covariance matrices of x t and n t . As the background is spatio-temporally constant, the covariance matrix of b is zero. Recalling that the emitters are uncorrelated by assumption, we deduce that R x is diagonal. We thus set r x := diag(R x ). Furthermore, by the i.i.d. assumption on n t , we have that R n = sI M 2 , where s ∈ R + and I M 2 is the identity matrix in R M 2 ×M 2

. Equation (4) can thus be re-written as:

r y = (Ψ Ψ)r x + sv I
where denotes the column-wise Kronecker product, r y is the column-wise vectorization of R y and v I = vec(I M 2 ). In order to estimate r x and the variance s and promote sparsity on r x , we introduce the non-convex 0 -type CEL0 regularizing penalty proposed in [START_REF] Soubies | A continuous exact 0 penalty (CEL0) for least squares regularized problem[END_REF] and consider the following minimization problem:

arg min rx≥0,s≥0 1 2 r y -(Ψ Ψ)r x -sv I 2 2 + Φ CEL0 (r x ; λ), (5)
where λ > 0 a positive regularisation parameter and the CEL0 penalty Φ CEL0 (•; λ) is defined by [START_REF] Soubies | A continuous exact 0 penalty (CEL0) for least squares regularized problem[END_REF]

Φ CEL0 (r x ; λ) = L 2 i=1 λ-ai 2 2 |(r x ) i | - √ 2λ ai 2 1 {|(rx)i|≤ √ 2λ a i } ,
with a i = (Ψ Ψ) i being the i-th column of Ψ Ψ.

The functional in ( 5) is continuous and non-convex, but has the same minimizers as the corresponding 2 -0 problem. Compared to SPARCOM [START_REF] Solomon | SPARCOM: Sparsity based super-resolution correlation microscopy[END_REF], solutions of (5) are sparser. Moreover, the estimation of s is expected to improve upon the quality of the estimation.

Intensity and background estimation

Solving (5) provides an estimation of r x from which we deduce the support of x denoted by Ω x = {i :

x i = 0} ⊂ 1, . . . , L 2 as Ω rx = Ω x .
We then estimate the intensity of x only on its support, and at the same time the spatially constant background b = b1 M 2 , b ≥ 0, by solving arg min

x∈R |Ωx| + , b≥0 1 2 y -Ψ Ω x -b1 M 2 2 2 + µ ∇ Ω x 2 2 ( 6
)
where the data term models the presence of Gaussian noise, µ > 0 is a regularization parameter, Ψ Ω ∈ R M 2 ×|Ωx| is a matrix whose i-th column is extracted from Ψ for index i ∈ Ω x and the regularization term is the squared norm of the discrete gradient restricted to Ω x , i.e.:

∇ Ω x 2 2 := 1 2 i∈Ωx j∈N (i)∩Ωx (x i -x j ) 2 ,
where N (i) denotes the 8-pixel neighbourhood of i ∈ Ω x .

ALGORITHMIC IMPLEMENTATION

We use alternating minimization to solve (5) (see Algorithm 1). For the estimation of r x , we follow [START_REF] Gazagnes | High density molecule localization for super-resolution microscopy using CEL0 based sparse approximation[END_REF] and use the iteratively reweighted 1 algorithm (IRL1). A good initialization r x 0 is the 1 -regularized solution. An explicit expression for s can be obtained from the (unconstrained) optimality condition. This estimate is then projected onto the set of positive solutions in a standard way. For solving [START_REF] Solomon | SPARCOM: Sparsity based super-resolution correlation microscopy[END_REF] we use again alternating minimization (see Algorithm 2) and solve each subproblem by a standard quadratic programming (QP) algorithm. Here only an initialization for the constant background value is required: a good choice for it is the median of y. For both algorithms, we consider stopping criteria based on the relative difference between consecutive iterates and on a maximum number of iterations.

Algorithm 1 COL0RME: Support estimation -1st step

Require: ry ∈ R M 4 , rx 0 ∈ R L 2 , s 0 ∈ R+, λ > 0 repeat compute weights ω r k x i ∈ ∂ΦCEL0(r k x ; λ) rx k+1 = arg min rx∈R L 2 + 1 2 ry -(Ψ Ψ)rx -s k v I 2 2 +λ L 2 i=1 ω i r k x |(rx)i| s k+1 = arg min s∈R + 1 2 ry -(Ψ Ψ)r k+1 x -sv I 2 2
until convergence return Ωx, s Algorithm 2 COL0RME: Intensity estimation -2nd step

Require: yt ∈ R M 2 , Ωx, b 0 , µ > 0 repeat x k+1 = arg min x∈R |Ωx| + 1 2 y -Ψ Ω -b k 1 M 2 2 2 + µ ∇Ωx 2 2 b k+1 = arg min b∈R + 1 2 y -Ψ Ω x k+1 -b1 M 2 2 2 until convergence return x, b

RESULTS AND DISCUSSION

Simulated data

We start by applying COL0RME to images of tubular structures simulating standard microscope acquisitions with standard fluorescent dyes. The spatial pattern is taken from the MT0 microtubules training dataset uploaded for the SMLM 2016 1 , see Fig. 1b. Intensities are obtained by using the SOFI simulation tool [START_REF] Girsault | SOFI simulation tool: A software package for simulating and testing super-resolution optical fluctuation imaging[END_REF]. Namely, we simulate temporal fluctuations and create videos of T = 100 and T = 700 frames at a frame rate of 100 frames per second (fps). The 1 http://bigwww.epfl.ch/smlm/datasets/index.html simulation parameters are set as follows: 20ms for on-state average lifetime, 40ms for off-state average lifetime and 20s for average time untill bleaching. The emitter density is equal to 10.7 emitters/pixel/frame, while the FWHM of the PSF is approximately 229nm. We create two different noisy datasets in order to evaluate the results of COL0RME and to compare them with the ones obtained by SRRF [START_REF] Gustafsson | Fast live-cell conventional fluorophore nanoscopy with ImageJ through super-resolution radial fluctuations[END_REF] and SPARCOM [START_REF] Solomon | SPARCOM: Sparsity based super-resolution correlation microscopy[END_REF], which, similarly, also exploit the temporal fluctuations of molecules. For the first dataset, we generate on average 1000 photons/frame per emitting molecule and b * = 100 photons/frame per pixel to simulate the out-offocus molecules, which we consider as background (BG). For the second dataset (which is more realistic), we set these values to 500 and b * = 2500, respectively. For both datasets we add Gaussian noise of 20dB to simulate the presence of electronic noise.

As we know the ground truth positions of the emitters, we can evaluate the localization precision of each method in terms of the Jaccard index, a quantity in the range [0, 1] computed as the ratio between correct detections (up to some tolerance δ > 0) and the sum of correct detections, false positives and false negatives. Table 1 reports the average Jaccard index (δ = 40nm) for 20 noise realizations. Note that, in this table, for both SRRF and SPARCOM the reported results have been obtained after a thresholding step to eliminate the many false positives in the reconstructed support, while COL0RME does not require any post-processing and achieves precise localization even in hard experimental conditions. In Fig. 1 we report the final reconstructions before thresholding showing that COL0RME better preserves the real structure of the tubulin. This is a different feature from SRRF whose reconstruction is overall visually more pleasant, but which creates artifacts as it cannot separate close-by filaments. COL0RME exhibits a more accurate reconstruction than SPARCOM whose result has many false positives around the tubulin. Thanks to the intensity estimation step, COL0RME is also able to estimate the brightness values which correspond to real signal intensities. The peak signal-to-noise ratio (PSNR) values w.r.t. to the ground truth intensity image and the estimated BG value b are provided in Table 2. We remark that the SRRF and SPARCOM do not actually estimate real signal intensities. Consequently, we had to adapt the range of their reconstructions according to the range of the ground truth image for the line profile graph in Fig. 1f.

Experimental data

COL0RME can also be applied for high-density acquisitions obtained by SMLM techniques. Even though in SMLM the molecules do not have a blinking behaviour, but rather a onto-off transition, we can consider as blinking the temporal behaviour of one pixel in high density videos with many molecules per pixel. We compared the methods COL0RME and SRRF on a patch extracted from a real dataset from the SMLM challenge 2013. frames and a FWHM of 351.8 nm. In Fig. 2 we observe that SRRF preserves better the broad structure of the specimen, while COL0RME better separates very close microtubules and does not create any background artifacts.

CONCLUSIONS

We proposed a super-resolution method named COL0RME, well suited for live-cell imaging and which can be easily applied to images obtained by common microscopes and conventional dyes. COL0RME takes advantage of the independent stochastic fluctuations of the fluorescent molecules and by solving a non-convex optimization problem in the covariance domain localizes the fluorophores with high precision. Differently from other super-resolution methods (SRRF, SPARCOM), COL0RME also includes an intensity estimation step, which is a valuable piece of information for biological analysis. We showed that COL0RME outperforms competing methods in terms of localization precision on both simulated and real data, and further computes an estimate of noise statistics and background information.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Results for the simulated 'High Background' dataset, with T = 700 before PP: (a) temporal average of the acquired stack (4× zoom), (b) ground truth, (c) COL0RME, (d) SRRF, (e) SPARCOM, (f) intensity profiles (SRRF and SPARCOM with adapted range) a b c

6 .

 6 COMPLIANCE WITH ETHICAL STANDARDS This work was conducted using biological data available in open access by EPFL SMLM datasets. Ethical approval was not required as confirmed by the license attached with the open access data.

Table 1 :

 1 The acquired video has T = 500 Jaccard index values with tolerance δ = 40nm for the different datasets. The results obtained by SRRF and SPARCOM are subject to a post-processing (PP) step to remove false positives.

				Low BG	High BG
		Method	T	100 700 100 700
		COL0RME (No PP) 0.51 0.66 0.28 0.51
		SRRF (PP)		0.5	0.5 0.47 0.47
		SPARCOM (PP)	0.32 0.39 0.17 0.31
	T	100	700		100	700
	PSNR 26.6	29.21		23.31	28.37
	b	92.3	102		2435	2454

Low BG, b * = 100 High BG, b * = 2500

Table 2 :

 2 COL0RME PSNR values (in dB) and estimated back-

ground b for the four datasets (average over 20 noise realizations).

ACKNOWLEDGMENTS

The work of VS and LBF has been supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002. LC acknowledges the support received by UCA IDEX JEDI and by the NoMADS RISE H2020 project 777826.