Probing the 3D Alpine lithosphere by Full Waveform Inversion of the AlpArray teleseismic data – LiSAlps

Full Waveform Inversion (FWI) has emerged since one decade as the baseline imaging method in many domains (exploration geophysics, earthquake seismology, medical imaging, nondestructive control). The overall goal is to estimate the constitutive properties of a medium from indirect measurements provided by elastic waves triggered by a layout of sources and recorded by a network of sensors. The mathematical model that describes propagation of elastic waves is the partial differential wave equation. The properties of the medium to image are gathered in the coefficients of the equation and the measurements over time collected during an experiment are the solution of this equation at the positions of the receivers for multiple right-hand sides (the sources). Therefore, FWI aims to estimate the spatially-varying coefficients of a partial differential equation from the parsimonious measurements of its solution. This is tackled by solving an inverse problem with local optimization techniques (gradient-based methods) by minimizing a function measuring a distance between the recorded data and the simulated counterpart (generally, the least-squares norm of the difference). The motivation behind FWI is to use the full information content of the data by exploiting each sample discretizing the seismograms (recording over time of waves triggered by a source and reaching a receiver) collected during the experiment. The add-value provided by this all-at-once approach is a theoretical resolution of the order of the wavelength and the ability to image all of the constitutive properties governing wave propagation. The price to pay is the nonlinearity of the inverse problem resulting from the oscillatory nature of waveforms, the ill-posedness of the imaging resulting from incomplete illumination of the Earth’s interior from the surface and the computational cost of full-wave imaging techniques. In teleseismic imaging, the sources (the earhquakes) are distant from the lithospheric target located below the network of stations. This setting requires specific modeling strategies where a first simulation in a simplified global Earth is performed to propagate a reference background wavefield from the source to the target before its injection at the boundaries and its propagation in the target toward the stations. The challenge of the inverse problem is to reach the desired resolution during imaging by extracting the information contained in the incident forward-scattered wavefield but also in the second-order back-scattered wavefield from the lithospheric heterogeneities after a first reflection from the topography. This raises multiple methodological challenges related to numerical simulation of wave propagation in large computational domains and robust optimization algorithms equiped with regularization able to manage incomplete and noisy data and reconstruct multiple classes of parameters with contrasted signatures in the data.

Before running FWI, it is crucial to perform a careful quality control of the data to prevent injecting outliers in the inversion. These outliers can have been generated by the failure of some stations or by an unfavorable radiation pattern of the source according to the azimuth angle of the event at the stations leading to a poor signal to noise ratio of the P wave at the station. A first dataset was selected gathering 20 events and the 600 stations of the AlpArray experiment according several criteria based on the phase, amplitude and polarization of the arrivals recorded by the multi-component sttaions. This dataset is based on the P-wave first arrival and its coda in a time window the length of which should be defined. The bandwidth used for inversion ranges between 10s and 30s with a limited window length of 15s. Wave simulation in the target is performed with the spectral element method on h-adaptive hexahedral mesh while the computatation of the reference wavefields in the global Earth is performed with AxiSEM software. The mesh of the targeted domain contains contains 126 900 elements parametrized with P4 Lagrange polynomials. The initial model for FWI is a smoothed version of the AK135 model including the topography of the Alps. The inversion is implemented with the quasi-Newton l-BFGS algorithm and 65 iterations are performed. The initial source signatures are extracted from the GlobalCMT catalogue and are refined on the fly during FWI in alternance with the model update, which is parametrized by the P and S wave velocities and density. Depth slices at 10km depth of the P-wave perturbation model allow us to delineate well-documented structures in the upper crust such as the large sedimentary basins (Po plain, Molassic basin, Rhone-Bresse basin) and the Ivrea body. Moreover, a first map of the Moho has been built, which reveals the crustal thickening beneath the Alps and the Apennines, a mild crustal thinning below the Rhin graben and the Rhone valley and a very thin crust of the Ligurian basin. Concerning upper mantle imaging, vertical sections of the P-wave model perturbation reveal the Alpine and Apennines slabs, whose track can be followed from 75km depth to 350km. This work is ongoing to produce the next generation of elastic models of the Alps by involving more events in the inversion, by increasing the window length, and by widening the period band.

The short-term pespectives of this work are to involve more events in the data, increase the recording length during inversion and push FWI at higher frequencies (0.3Hz) to improve the FWI resolution while checking the acquisition footprint in the reconstructed models (a limiting factor for pushing FWI at high frequencies is the detrimental effect of sparse acquisitions in terms of aliasing). This will conclude the first phase of the project with a new generation elastic isotropic model of the crust and upper mantle of the Alps. In parallel, a S-wave dataset will be built to improve the results inferred from the P wave and its coda; in particular we expect significant improvement of the Vs model when the incident S waves are involved in the inversion. This will be also the first step toward anisotropic FWI where we will be able to account for shear-wave splitting. On the application side, the second phase will focus on the Ligurian knot where the data of the cifalps2 experiment will supplement those of AlpArray. The same numerical strategies as those used durant the first phase will be used except that the FWI will be pushed up to a frequency of 1Hz. From the methodological viewpoint, we will assess non conventional regularization techniques allowing us to conciliate a high-resolution imaging of contrasted structures in the crust with the imaging of smoother structures in the upper mantle. This can be achieved with hybrid regularizations combining total variation regularization suitable for blocky media and Tikhonov regularization suitable for smooth media. These regularizations should be combined mathematically such that the regularization can be tailored to the local statistical properties of the medium. The l1-norm based regularization are also useful to mitigate the footprint of sparse acquisition on imaging. To optimally tune the regularization for this purpose, it will be useful to build a geomodel of the Alps to test the regularized FWI on a fully controled experiment (in terms of acquisition, data and model). Finally, new formulations of FWI were recently proposed where the accuracy of the wavefields computed in the background model is improved by adding a feedback term to the data in the source. We call these wavefields data-assimilated wavefields. These wavefields allow us to build more accurate sensitivity kernel of FWI during parameter updating. Although this method has been more specifically designed to deal with strong nonlinearity (which is not the central issue in earthquake seismology due to the availability of low frequencies), this should not prevent to assess the added-value of this new FWI formulation on teleseismic applications.

Ongoing

The Alpine chain has historically been one of the main natural laboratories to study orogenesis due to its outstanding outcrop conditions. Despite countless investigations, the deep structures imaged by geophysical methods and related tectonic processes remain controversial, due to the high complexity of the chain but also from the lack of high-resolution structural information at depth. To fill this gap, the European AlpArray consortium deployed 628 broadband stations, spaced less than 52 km apart, across the whole chain, and 30 ocean-bottom seismometers in the Ligurian basin, hence providing a unique opportunity for a step change in the 3D imaging of the Alpine lithosphere and asthenosphere. The LisAlps project proposes to apply Full Waveform Inversion (FWI) on the teleseismic data recorded during AlpArray to build (1) a new reference high-resolution multi-parametric (1500x700 km) model of the alpine lithosphere and asthenosphere down to 700km depth from the entire network and a catalogue of ~300 teleseismic earthquakes (periods: 5s-20s); (2) a high-resolution model of the lithosphere in the western Alps around the structurally-complex Ligurian knot. Through this imaging, we want to: [a] decipher the debated geometry and petrology of the continental subduction between the European and Adriatic plates, [b] sharpen the location of the seismicity in the western Alps, in particular the deep events beneath the Po plain, [c] extend at greater depths with an increased resolution the imaging of the continental subduction to test the lateral and downdip continuity of the slabs between the western and central Alps and refine the slab structure in the transition zone between the Alps and the Apennines, [d] image the interaction between opposite verging Alpine and Apennine slabs and assess the possible role of such interaction in the compression within Northern Apennines, [e] understand the compressional reactivation and ensuing Oligo-Miocene inversion of the northern Ligurian margin, [f] Solve the origin of the high-rate microseismicity and occasional strong earthquakes (e.g., 1887 Mw6.9 event) over the North Ligurian domain for geohazard assessment in the coastal areas of the French Riviera and Liguria. Answering these questions raise at least three methodological challenges. Unlike other methods used in seismology, FWI exploits both the amplitude and phase of all the arrivals recorded in a continuous time window. This ability provides a wavelength-scale resolution and the sensitivity to several classes of physical properties such as P and S wavespeeds, density, attenuation and anisotropy. However, the coarse sampling of the sources and stations in earthquake seismology can generate aliasing artefacts. To mitigate them, we propose to assess compress sensing via sparsity-promoting regularization. A second challenge is the estimation of “second-order” parameters as anisotropy and attenuation, which are two important proxy for several Earth's properties such as temperature, composition (serpentinization, melt), state of stress and asthenospheric mantle flows. A third challenge is the non linearity of the inversion induced by fitting oscillating signals. To tackle this issue, we propose to assess a new “upside down” FWI paradigm, which reconstructs “data-assimilated” wavefields that fit the observations through a relaxation of the wave equation before estimating the Earth's parameters from these wavefields by minimizing the wave equation errors (a linear problem). The better data-assimilated wavefields mimic the true unknown wavefields, the more the parameter estimation is fast and robust. We believe that the full down-top/top-down illumination of the lithosphere provided by the specific teleseismic setup should foster accurate data-assimilated wavefield reconstruction. The proof of concept of the approaches we propose has shown promising results. We propose to assess them against the real case studies scheduled during LisAlps.