\makeatletter
\@ifundefined{HCode}
{\documentclass[screen, CRGEOS, Unicode, Thematic,published]{cedram}
%\usepackage[sort&compress,square,comma,authoryear]{natbib}
\usepackage[square,comma,authoryear]{natbib}
\newenvironment{Table}{\begin{table}}{\end{table}}
\renewcommand\bibsection{\section*{\refname}}
\newenvironment{noXML}{}{}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}}
\def\tabnote#1{\vskip4pt\parbox{.92\linewidth}{#1}}
\def\xxtabnote#1{\vskip4pt\parbox{.94\linewidth}{#1}}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\botline{\\\hline}
\usepackage[figuresright]{rotating}
\RequirePackage{etoolbox}
\usepackage{upgreek}
\def\jobid{crgeos20230948}
%\graphicspath{{/tmp/\jobid_figs/web/}}
\graphicspath{{./figures/}}
\newcounter{runlevel}
\usepackage{multirow} 
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\let\MakeYrStrItalic\relax
\def\ndash{\text{--}}
\def\refinput#1{}
\def\back#1{}
\def\ubreak{\break}
\csdef{Seqnsplit}{\\}
\DOI{10.5802/crgeos.261}
\datereceived{2023-10-25}
\daterevised{2024-04-19}
\dateaccepted{2024-04-26}
\ItHasTeXPublished
}
{
\PassOptionsToPackage{authoryear}{natbib}
\documentclass[crgeos]{article}
\usepackage[T1]{fontenc} 
\def\CDRdoi{10.5802/crgeos.261}
\def\xmorerows#1#2{\morerows{#1}{#2}}
\let\newline\break
\def\xxtabnote#1{\tabnote{#1}}
\let\ubreak\relax
\makeatletter
\def\CDRsupplementaryTwotypes#1#2{}
\def\selectlanguage#1{} 
}
\makeatother

\usepackage{float}

\dateposted{2024-09-25}
\begin{document}

\begin{noXML}

\CDRsetmeta{articletype}{research-article}

%\makeatletter
%\def\TITREspecial{\relax}
%\def\cdr@specialtitle@english{New Developments in Passive Seismic Imaging and Monitoring}
%\def\cdr@specialtitle@french{Nouveaux d\'eveloppements en mati\`ere d'imagerie sismique passive et de surveillance}
%\makeatother

\title{Methodological advances in seismic noise imaging of the Alpine area}

\alttitle{{Avanc\'{e}es m\'{e}thodologiques sur l'imagerie du bruit
sismique dans la zone alpine}}

\author{\firstname{Anne} \lastname{Paul}\CDRorcid{0000-0002-5973-1786}\IsCorresp}
\address{Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, UGE,
ISTerre, Grenoble, France}
\email[A. Paul]{anne.paul@univ-grenoble-alpes.fr}

\author{\firstname{Helle A.} \lastname{Pedersen}\CDRorcid{0000-0003-2936-8047}}
\addressSameAs{1}{Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS,
IRD, UGE, ISTerre, Grenoble, France}
\email[H. A. Pedersen]{helle.pedersen@univ-grenoble-alpes.fr}

\author{\firstname{Thomas} \lastname{Bodin}\CDRorcid{0000-0002-7393-6112}}
\address{Univ. Lyon, Univ. Lyon 1, ENSL, UJM-Saint-Etienne, CNRS,
LGL-TPE, F-69622, Villeurbanne, France}
\email[T. Bodin]{thomas.bodin@ens-lyon.fr}

\author{\firstname{Emanuel} \lastname{K\"{a}stle}\CDRorcid{0000-0002-1866-2323}}
\address{Institute of Geological Sciences, Freie Universit\"{a}t
Berlin, Berlin, Germany}
\email[E. K\"{a}stle]{emanuel.kaestle@fu-berlin.de}

\author{\firstname{Dorian}~\lastname{Soergel}\CDRorcid{0000-0003-0847-1635}}
\addressSameAs{1}{Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS,
IRD, UGE, ISTerre, Grenoble, France}
\curraddr[D. Soergel]{Berkeley Seismological Laboratory, University of
California, Berkeley, CA~94720, USA}
\email[D. Soergel]{soergeldorian@gmail.com}

\author{\firstname{Chlo\'{e}} \lastname{Alder}\CDRorcid{0000-0001-6093-2772}}
\addressSameAs{2}{Univ. Lyon, Univ. Lyon 1, ENSL, UJM-Saint-Etienne,
CNRS, LGL-TPE, F-69622, Villeurbanne, France}
\email[C. Alder]{chloe.alder@gmx.fr}

\author{\firstname{Yang} \lastname{Lu}\CDRorcid{0000-0002-5192-8840}}
\addressSameAs{1}{Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS,
IRD, UGE, ISTerre, Grenoble, France}
\curraddr[Y. Lu]{Department of Meteorology and Geophysics, University
of Vienna, Vienna, Austria}
\email[Y. Lu]{yang.lu@univie.ac.at}

\author{\firstname{Ahmed} \lastname{Nouibat}\CDRorcid{0000-0002-3069-2701}}
\addressSameAs{1}{Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS,
IRD, UGE, ISTerre, Grenoble, France}
\curraddr[A. Nouibat]{ITES-EOST, CNRS/Universit\'{e} de Strasbourg, 67000
Strasbourg, France}
\email[A. Nouibat]{ahmed.nouibat@univ-grenoble-alpes.fr}

\shortrunauthors

\CDRGrant[ANR]{ANR-15-CE31-0015}
\CDRGrant[ANR]{ANR-10-LABX-56}
\CDRGrant[EU]{716542}
\CDRGrant[ANR]{ANR-20-CE49-0007}
\CDRGrant[DFG]{2403/21-1}

\keywords{\kwd{Seismic tomography}\kwd{Ambient noise}\kwd{Bayesian
inversion}\kwd{Seismic anisotropy}\kwd{Lithospheric
structure}\kwd{European Alps}}

\altkeywords{\kwd{Tomographie sismique}\kwd{Bruit ambiant}
\kwd{Inversion Bay\'esienne}\kwd{Anisotropie sismique}
\kwd{Structure lithosph\'erique}\kwd{Alpes}}

\thanks{Agence Nationale de la Recherche, France (contract
ANR-15-CE31-0015), Labex OSUG@2020 (Investissement d'Avenir,
ANR-10-LABX-56), European Research Council under the European Union
Horizon 2020 research and innovation program (grant agreement no.
716542 -- TRANSCALE), project LisAlps (contract ANR-20-CE49-0007),
German Science Foundation DFG (SPP-2017, Project Ha 2403/21-1).}

\begin{abstract}
Methodological advances in seismic tomography are often driven by the
quality of data sets. The dense and homogeneous spatial coverage of the
AlpArray seismic network, including hundreds of permanent and temporary
broadband stations, has motivated a series of methodological
developments of ambient-noise-based tomography of the lithosphere
across the entire Alps-Apennines regions, which have been published and
are reviewed here. To take full advantage of the ocean-bottom
seismometers (OBS) in the Ligurian-Provence basin, reconstructed
Rayleigh wave signals between OBS have been improved by second-order
correlations with onland stations. A Bayesian or fully transdimensional
formalism has been introduced in both steps of isotropic ambient noise
tomography. The three-dimensional S-wave velocity models have been
further improved by wave-equation based inversions accounting for the
physics of seismic wave propagation, including elastic--acoustic
coupling at the sea bottom. A beamforming approach has been developed
to avoid systematic errors in the measurement of azimuthal anisotropy
from seismic noise. Probabilistic inversions for depth variations of
azimuthal and radial anisotropy have provided robust estimates of
anisotropic parameters in the crust and upper mantle that differ
significantly from earlier surface-wave tomography studies.
These methodological improvements have taken the full benefit of the
quality of available seismic data to significantly improve knowledge of
the seismic structure of the crust and shallow mantle beneath the
Alps-Apennines system. Our findings include detailed mapping of strong
and abrupt Moho depth changes under the Western Alps, contrasting
orientations of fast velocity directions between the upper and lower
Alpine crust, and the absence of significant radial anisotropy
everywhere in the European crust and shallow upper mantle, except in
the Apenninic lower crust. These methods can be applied to similar
dense arrays with equivalent potential benefits.
\end{abstract}


\begin{altabstract}
Les progr\`{e}s m\'{e}thodologiques dans le domaine de la tomographie
sismique d\'{e}pendent souvent de la qualit\'{e} des ensembles de
donn\'{e}es. La couverture spatiale dense et homog\`{e}ne du r\'{e}seau
sismique AlpArray, comprenant des centaines de stations permanentes et
temporaires \`{a} large bande, a motiv\'{e} une s\'{e}rie de
d\'{e}veloppements m\'{e}thodologiques de la tomographie de la
lithosph\`{e}re bas\'{e}e sur le bruit ambiant dans l'ensemble des
r\'{e}gions Alpes-Apennins, qui ont \'{e}t\'{e} publi\'{e}s et sont
pass\'{e}s en revue ici. Pour tirer pleinement parti des
sismom\`{e}tres de fond de mer (OBS) dans le bassin
Liguro-Proven\c{c}al, les signaux d'ondes de Rayleigh reconstruits
entre les OBS ont \'{e}t\'{e} am\'{e}lior\'{e}s par des
corr\'{e}lations de second ordre avec les stations terrestres. Un
formalisme bay\'{e}sien ou enti\`{e}rement transdimensionnel a
\'{e}t\'{e} introduit dans les deux \'{e}tapes de la tomographie
isotrope du bruit ambiant. Les mod\`{e}les tridimensionnels de vitesse
des ondes S ont \'{e}t\'{e} am\'{e}lior\'{e}s par des inversions
bas\'{e}es sur l'\'{e}quation des ondes qui tiennent compte de la
physique de la propagation des ondes sismiques, y compris le couplage
\'{e}lastique-acoustique en fond de mer. Une approche par formation de
voie (beamforming) a \'{e}t\'{e} d\'{e}velopp\'{e}e pour \'{e}viter les
erreurs syst\'{e}matiques dans la mesure de l'anisotropie azimutale
\`{a} partir du bruit sismique. Des inversions probabilistes pour
obtenir les variations en profondeur de l'anisotropie azimutale et
radiale ont fourni des estimations robustes des param\`{e}tres
anisotropes dans la cro\^{u}te et le manteau sup\'{e}rieur qui
diff\`{e}rent sensiblement des \'{e}tudes ant\'{e}rieures de
tomographie par ondes de surface. Ces am\'{e}liorations
m\'{e}thodologiques ont permis de tirer pleinement parti de la
qualit\'{e} des donn\'{e}es sismiques disponibles pour am\'{e}liorer de
mani\`{e}re significative la connaissance de la structure sismique de
la cro\^{u}te et du manteau peu profond sous le syst\`{e}me
Alpes-Apennins. Nos r\'{e}sultats comprennent une cartographie
d\'{e}taill\'{e}e des changements de profondeur importants et abrupts
du Moho sous les Alpes occidentales, des orientations contrast\'{e}es
des axes de vitesse rapide entre la cro\^{u}te alpine sup\'{e}rieure et
inf\'{e}rieure, et l'absence d'anisotropie radiale significative
partout dans la cro\^{u}te europ\'{e}enne et le manteau sup\'{e}rieur
peu profond, \`{a} l'exception de la cro\^{u}te inf\'{e}rieure des
Apennins. Ces m\'{e}thodes peuvent \^{e}tre appliqu\'{e}es \`{a} des
r\'{e}seaux denses similaires avec des avantages potentiels
\'{e}quivalents.
\end{altabstract}

\maketitle

\vspace*{6pt}

\twocolumngrid

\end{noXML}

\section{Introduction}\label{sec1}

The Alpine mountain range is part of the continental collision ranges
created by the convergence of the Eurasian and African plates in the
Mediterranean region. It results from the subduction of the Alpine
Tethys under the Adriatic microplate since the Late Cretaceous, and the
subsequent continental collision between the European and Adriatic
paleomargins in the Cenozoic \citep[e.g.,][]{Handyetal2010}.
The Alpine belt is an outstanding example of subduction and continental
collision studies, which has been investigated by geologists for more
than 150 years. Several major concepts of modern geology have been
developed in the Alps, such as nappes, when the prominence of
horizontal over vertical displacements was proposed by Emile Argand
\citep{Argand1922}. More recently, the identification of coesite, which
is a high-pressure polymorph of quartz, in gneisses of the Dora Maira
massif (south-western Alps, Italy) led \citet{Chopin1984} to propose that
continental crust may be subducted to depths of 90 km or more. The
amount of geological knowledge about the Alps is unparalleled in any
other mountain range, and they provide a unique natural laboratory to
advance our understanding of orogenesis and its relationship to present
and past mantle dynamics. The Alpine mountain belt is also a populated
area where millions of Europeans are affected by its topography,
geology and associated natural hazards such as earthquakes or
landslides. Yet, accurate information on the lithospheric structure of
that emblematic and populated mountain range was hampered by
insufficient and spatially heterogeneous geophysical data until the
last few years. Filling that gap was the primary motivation for the
AlpArray initiative, which gathered a large number of European research
institutions to deploy a dense and homogeneous temporary broadband
seismological network over the Alps and its forelands to complement the
permanent networks~\citep{Hetenyietal2018a}.\looseness=-1

In addition to gaining knowledge of Alpine lithospheric structure, the
high spatial coverage of AlpArray temporary stations and permanent
networks has given us the opportunity to develop new methods of seismic
tomography at this large regional scale based on ambient noise. A
review of these methodological developments is gathered in this paper
because they could be applied with great benefit to other similary dense
networks.

\subsection{Seismic imaging in the Alps}

Seismic imaging is an essential complement to geological studies to
build lithospheric-scale interpretive models and improve the
understanding of the dynamics of the mountain belt in space and time.
The geometry and depth of the crust-mantle boundary, only accessible
with geophysics and in particular active-source seismology and
earthquake seismology, are key information for geological and geodynamic modeling.
Each of these methods has limitations and provides partial information,
which is why noise-based tomography methods are a valuable complement.
In the Alps, deep seismic sounding (DSS) experiments including
ECORS-CROP in the Western Alps \citep[e.g.,][]{Nicolasetal1990}, NFP-20 in
the Central Alps \citep[e.g.,][]{Freietal1990}, and TRANSALP in the
Eastern Alps \citep[e.g.,][]{Luschenetal2006} provided crucial data for
interpretive crustal-scale sections along a number of crooked lines
mostly transverse to the belt. The results from these localised studies
cannot be extrapolated to other locations along the belt due to its
arcuate, non-cylindrical geometry. DSS profiles provide high-frequency
reflectivity images of the crust, hence sharp and sub-horizontal
velocity contrasts with no or poor information on absolute velocities. 
Thus, they cannot be interpreted in terms of \mbox{petrology.} Moreover, the
European Moho, defined as the base of the highly reflective lower crust,
was not detected below the intensely deformed and highly heterogeneous
crust of the internal zones in the ECORS-CROP and NFP-20 W5 and E1
profiles across the Western and Central Alps
\citep[e.g.,][]{Nicolasetal1990,MarchantStampfli1997}. The deep European
Moho has only been detected by the ECORS-CROP wide-angle seismic
reflection experiment under part of the internal zones to a maximum
depth of 55~km \citep{ECORS-CROPDeepSeismicSoundingGroup1989}.
\looseness=1

A similar type of information on velocity contrasts beneath seismic
stations is provided at lower frequencies by receiver function (RF)
analysis in particular for Moho depth estimates 
\citep[e.g.,][]{Kummerowetal2004,Lombardietal2008,Zhaoetal2015,Hetenyietal2018b,Pauletal2022,Michailosetal2023}. 
In some specific
cases, such as in the presence of a very heterogeneous, scattering
crust, receiver function analysis can be more effective at detecting
the crust-mantle boundary than reflection seismics because it uses
low-frequency, high-energy waves from teleseismic earthquakes that
travel across the heterogeneous crust only once. Moreover, receiver
functions computed for station arrays provide estimates of Moho depth
with 2-D coverage. Indeed, \citet{Spadaetal2013} generated a Moho
depth map of Italy, including the Alpine region, which combines DSS
data with receiver functions. \citet{Spadaetal2013}'s model
displays three Moho surfaces, European, Adriatic-Ionian and
Ligurian--Corsican--Sardinian--Tyrrhenian. As this model allows only one
Moho at a given location, the Moho surfaces never overlap, even though
the European lithosphere is known to underthrust Adria.


The polarity of converted waves in receiver functions provides useful
information on the sign of velocity change with depth. For example,
\citet{Zhaoetal2015} observed negative-polarity P-to-S converted waves
in RF of the CIFALPS profile close to the so-called Ivrea body positive
Bouguer anomaly (blue line in Figure~\ref{fig1}). This gravity anomaly
high is known since the first geophysical experiments in the Alps that
also reported high-velocity refracted waves ($V_{\mathrm{p}}=7.4$~km/s) at
10~km depth in the same area of the Italian Piemonte region
\citep{ClossLabrouste1963}. The source of this gravity and seismic
velocity anomaly is called the Ivrea body and it is\break interpreted as a
slice of Adriatic upper mantle at unusually shallow depth
\citep[e.g.,][]{Nicolasetal1990}. The negative-polarity converted phases
in the RF of the CIFALPS profile were the first evidence for the
inverted Moho beneath the Ivrea body, as the contact between the
high-velocity Adriatic mantle wedge on top and the lower velocity
European crust or subduction interface below \citep{Zhaoetal2015}.
These negative-polarity converted waves from 20--60~km depth and the
positive-polarity conversions at 75--80~km depth were the first seismic
evidence for continental subduction of the European lithosphere beneath
Adria, according to \citet{Zhaoetal2015}. Like deep seismic sounding,
however, receiver functions yield no clues on absolute velocities.


\begin{figure*}
{\vspace*{-5pt}}
\includegraphics{fig01}
{\vspace*{-7pt}}
\caption{\label{fig1}Maps of broadband seismic stations in the Alps and
surrounding regions including all stations used in the studies
discussed in Sections~\ref{sec2}--\ref{sec3}; (a)~before 2010: red
triangles are public permanent stations in 2009; (b)~after 2020: red
triangles are permanent stations in 2019; green triangles are temporary
stations used in the studies described in this paper (EASI, AlpArray,
CIFALPS and CIFALPS-2). The blue line in (a) and (b) is the 0 mgal
contour of the Ivrea Bouguer anomaly. Black lines: main geological
units and faults from the tectonic map compiled by M.~R. Handy. The
thick black line CC$'$ is the CIFALPS profile used in Figure~\ref{fig3}.
DM: Dora Maira massif; IGA: Ivrea gravity anomaly.}
\end{figure*}


Local earthquake tomography (LET), which relies on the inversion of
body-wave arrival times (mostly direct P and S waves) from local
earthquakes for absolute velocities ($V_{\mathrm{p}}$ and $V_{\mathrm{s}}$ or $V_{\mathrm{p}}/V_{\mathrm{s}}$) is an
efficient tool to get 3-D images and complement 2-D DSS or receiver
function reflectivity images. The size of the imaged crustal volume and
the resolution of the tomography depend on the distribution of seismic
stations at the surface and earthquake hypocenters at depth. In the
Alps, the seismicity level is low to moderate, with rare events of
magnitude ${>}$5. Many years of recording are therefore required to reach
sufficient ray coverage in local earthquake tomography studies of the
whole Alpine region like in \citet{Diehletal2009}: 12 years
(1996--2007) for 1500 events of $M\textsubscript{l}> 2.5$. Another issue is that most
earthquake sources in the Alps are shallow (focal depths $<$ 15 km),
except in the westernmost part of the Po basin where hypocenters on the
Rivoli-Marene fault reach 70--80 km depth \citep{Evaetal2015}. Due to
low magnitudes, these deep events are only detected by stations close
to epicenters. They, however, proved useful to image the deep parts of
the subduction wedge of the western Alps by LET including the
high-velocity Ivrea body
\citep{Solarinoetal1997,Pauletal2001,Solarinoetal2018,Virieuxetal2024}.
Imaging the crust to Moho depth with LET in most of the Alps requires
records of Pn waves, thus sufficient-magnitude earthquakes and station
arrays of large spatial extent. Until recent years, this last condition
was only reached by compiling records of a large number of regional and
national networks with heterogeneous data sharing policies and
traveltime picking strategies, at the cost of heavy homogenization work
[13 networks used in \citealp{Diehletal2009,Bagaglietal2022}, and
\citealp{Virieuxetal2024}]. \looseness=-1

In addition to low-to-moderate seismicity and heterogeneous coverage of
permanent seismic networks, the strong crustal heterogeneity and its
very 3-D character are other challenges in seismic tomography of the
Alpine lithosphere. Crustal heterogeneity is not surprising for a
collisional belt with a long tectonic history. However, the arcuate
shape of the western Alps, the presence of large and deep basins
including the Po basin and its unconsolidated sediment layers and high
level of anthropogenic noise close to the heart of the belt make
seismic imaging of the Alpine lithosphere particularly challenging.

Since \citet{Akietal1977}, the most basic and often used method to
image seismic heterogeneities in the upper mantle is teleseismic
tomography. It is based on the inversion of observed relative arrival
times of P (or S) waves generated by earthquakes at teleseismic
distances (${>}$20\textdegree) for relative variations of P (or S) wave
velocity ($V_{\mathrm{p}}$, or $V_{\mathrm{s}}$) in the upper mantle beneath a seismic array.
Teleseismic tomography has been widely used to image fast-velocity
slabs interpreted as subducted continental and/or oceanic lithosphere
in the Alpine upper mantle
\citep[e.g.,][]{Lippitschetal2003,PiromalloMorelli2003,Zhaoetal2016a,Paffrathetal2021}. 
Contrasting results of these mantle tomography studies have led to
controversies about the geometry of the slabs at depth, including
whether the European slab is attached or detached in the western Alps
\citep[e.g.,][]{Zhaoetal2016a}, or whether Europe is the upper or lower
plate in the eastern Alps \citep[e.g.,][]{Lippitschetal2003}. These
questions are still of the utmost importance to understand the past and
present dynamics of the mountain belt. A major issue for teleseismic
tomography is the difficulty in separating the contribution of the
crust from that of the mantle in arrival times for near-vertical ray
paths \citep[e.g.,][]{Waldhauseretal2002}. This is of particular
importance in mountain ranges such as the Alps due to strong changes in
crustal thickness. 

The mantle structure can also be imaged with teleseismic surface-wave
tomography \citep[e.g.,][]{El-Sharkawyetal2020}. Frequency-dependent
traveltime data of surface waves are inverted for S-wave velocity as
for ambient-noise tomography, but using periods ${>}$30~s. To overcome the
poor sensitivity of such long period surface waves to crustal
structure, \citet{Kastleetal2018} have jointly inverted dispersion data
from\break ambient noise correlations in the short-period band (8--30~s) with
data from teleseismic surface-wave records at longer periods. The
horizontal resolution of the resulting shear-wave velocity models
however remains rather low in the upper mantle, limiting the value of
this type of model in the debates on slab geometry beneath the Alpine
range.

\subsection{Ambient noise imaging in the Alps before \mbox{AlpArray}}

Ambient noise tomography (ANT) is particularly well suited to imaging
the Alpine crust as a complement to the methods outlined in the
previous section because (a)~it does not require local earthquakes
\citep{Shapiroetal2005}, and (b)~the period range of seismic noise is
adequate for crustal imaging, in contrast to teleseismic surface wave
tomography which is dominated by longer wavelengths. A sufficient
coverage of the study region by seismic arrays is the only requirement
for a fairly resolved 3-D velocity model since each station becomes a
wave source for all other stations in the noise cross-correlation
process \citep{CampilloPaul2003,ShapiroCampillo2004}. 

Since the first application of ambient noise tomography to the Alps by
\citet{Stehlyetal2009}, station coverage has improved significantly in
density but also in spatial homogeneity. This improvement is due to
numerous new permanent broadband stations (in Austria, Germany, France,
etc.) and to temporary networks, most importantly the AlpArray seismic
network {[}AASN; \citealp{AlpArraySeismicNetwork2015}]. Indeed the
AASN, which included more than 600 stations, was designed to fill in
gaps between permanent stations and to homogenize spatial coverage in a
way that no location in the Alps was more than 30~km away from a
seismic station onland \citep{Hetenyietal2018a}. The AASN also had a
marine component with 29 ocean-bottom seismometers (OBS) deployed for 8
months in the Ligurian basin. It was complemented by denser
quasi-linear or 2-D temporary arrays on targets of specific interest.
The most important ones were EASI across the Eastern Alps
\citep{AlpArraySeismicNetwork2014}, CIFALPS and CIFALPS-2 across the
Western Alps \citep{Zhaoetal2016b,Zhaoetal2018}, and Swath-D in the
Central and Eastern Alps \citep{Heitetal2017}. Figure~\ref{fig1} shows the
seismological stations prior to 2010, and the stations which were used
in the studies\break summarised in this review. For seismic tomography in
general, and ambient noise tomography in particular, there is clearly a
time before and a time after\break AlpArray. 

Before AlpArray, the limited coverage of large parts of the Alpine area
only allowed for ambient noise tomography studies of the central Alps,
relying mainly on the dense Italian and Swiss permanent networks
\citep{Stehlyetal2009,Verbekeetal2012,Molinarietal2015}. The pioneering
application of ANT to the Alps by \citet{Stehlyetal2009} used records
of 150 broadband stations in Switzerland and neighbouring countries.
They applied a classical two-step procedure, with an inversion of
dispersion measurements for group velocity maps of Rayleigh waves in
the 5--80~s period band, followed by a non-linear Monte Carlo inversion
of the dispersion curve in each cell for a 1-D shear-wave velocity
model. The set of 1-D models were then merged into a so-called 3-D $V_{\mathrm{s}}$
model that should more properly be called pseudo 3-D. This study
provided a data constrained Moho depth map of the Western Alps, because
crustal thickness was a free parameter in the inversion. This Moho map
shared many similarities with the reference map computed by
\citet{Waldhauseretal1998} from depth-migrated controlled-source
seismic data, including strong and abrupt Moho depth changes from
25--30~km beneath the European forelands and the Po plain, to 55 km
beneath the internal Alpine arc from southern Switzerland to the Dora
Maira massif (Piemonte, Italy). This first successful application of
ANT in the Alps was considered a proof of concept. The tomography by
\citet{Stehlyetal2009} was followed by those of \citet{Verbekeetal2012}
and \citet{Molinarietal2015}, who expanded the station array to the
Italian and Slovenian permanent broadband networks and the tomography
to the Apennines, and used both group and phase velocity dispersion
data. 

In the last large-scale ANT before AlpArray, \citet{Kastleetal2018}
used records of 313 permanent stations covering the broad Alpine region
and the Apennines. To overcome the sparse station coverage in the
external, western and northern Alps, they complemented ambient noise
phase dispersion measurements with two-station measurements from
regional and\break teleseismic earthquake records. The ambient noise and
earthquake-based dispersion datasets agreed well enough in the
overlapping period band 8--60~s to be jointly inverted for Rayleigh and
Love wave phase velocity maps in the broad period range 4--250~s
\citep{Kastleetal2016}. In a second stage, \citet{Kastleetal2018}
jointly inverted Rayleigh and Love-wave phase dispersion data in each
cell for 1-D $V_{\mathrm{p}}$ and $V_{\mathrm{s}}$ models of the crust and upper mantle. According
to the authors, the 3-D $V_{\mathrm{s}}$ model of the upper mantle derived from joint
inversion has much higher resolution than when each individual dataset,
ambient noise and earthquake-based, is inverted alone. By averaging the
crustal thickness of the 500 best-fitting $V_{\mathrm{s}}$ models,
\citet{Kastleetal2018} computed a Moho depth map that compares
remarkably well with receiver function and DSS Moho depth estimates
along numerous profiles across the Alps and Apennines
\citep{Kummerowetal2004,Spadaetal2013,Zhaoetal2015}. 

The installation of the AlpArray temporary seismic network started in
Austria in the summer of 2015 and ${\sim}$200 of the planned 267
temporary stations were operating by mid-2016 \citep{Hetenyietal2018a}.
By adding the first months of AASN recordings (until June 2016) to a
four-year noise dataset from European-wide permanent broadband
networks, \citet{Luetal2018} achieved the first ANT of the broad Alpine
region with fairly homogeneous coverage and an average station spacing
of ${\sim}$50~km. This study was the first in a series of noise-based
isotropic and anisotropic tomography studies on the Alpine lithosphere
that have built on the spatial homogeneity and density of the AlpArray
dataset (from both permanent and temporary stations) to develop and
apply new data analysis and imaging methods. 

Indeed, the high spatial density of AASN and associated permanent and
temporary networks has driven methodological advances in ambient noise
tomography, just as the USArray Transportable Array has driven the
advent of Eikonal tomography \citep{Linetal2009}. The present paper
focuses on key methodological advances on isotropic
(Section~\ref{sec2}) and anisotropic (Section~\ref{sec3}) imaging. We
also show how the application of these new methods to Alpine data opens
new perspectives (Section~\ref{sec4}) for the geological and geodynamic
modeling of the Alpine belt.

\section{Isotropic ambient noise tomography}\label{sec2}

The homogeneous coverage provided by the AASN motivated a number of
ambient noise isotropic tomography studies on regional targets in and
around the Alps
\citep[e.g.,][]{Guerinetal2020,Molinarietal2020,Sadeghi-Bagherabadietal2021,Schippkusetal2018,Szanyietal2021}. 
In this section we focus on noise-based isotropic tomography studies of
the broad Alpine region and surroundings using all or most available
stations in Western Europe at the time of study. We describe three key
improvements to isotropic imaging with cross correlations in this
highly heterogeneous area. 

\subsection{Second-order cross-correlation techniques}

The Ligurian basin, which separates the Corsica--Sardinia block from the
southern coast of France and north-western Italy (Figure~\ref{fig1}) is
a back-arc basin generated by the rollback of the Adriatic slab in
Oligo-Miocene time \citep[e.g.,][]{Rolletetal2002}. Situated at the
transition between the Alps and Apennines mountain ranges, it was
considered an important target for the AlpArray temporary seismic
experiment, which therefore included 24 broadband ocean-bottom
seismometers (OBS, Figure~\ref{fig1}).

Ambient noise tomography of an offshore area is possible without OBSs
provided that it is surrounded by land stations. For example,
\citet{Magrinietal2022} have computed an S-wave velocity model of the
west-central Mediterranean, including the Ligurian basin, from
surface-wave tomography using records of onland stations in Europe and
North Africa. Lateral resolution in the Ligurian basin and its margins
remains limited, due to lack of OBS data, illustrating that the
integration of land and sea-based observations is a key target for
noise-based imaging in coastal and offshore areas. Such land-sea data
integration was the target of a specific effort for methodological
development by \citet{Nouibatetal2022b}. This work aimed at solving the
issues of ANT from OBS data, as OBSs are generally deployed for less
than a year, and signal quality is lower than that of land stations. In
this section, we summarize how \citet{Nouibatetal2022b} enhanced the
signal-to-noise ratio of noise correlations between sea-bottom stations
by computing second-order correlations with on-land stations as virtual
sources.

At periods longer than 20~s, OBS recordings are affected by tilt and
compliance noise induced by the soft seabed on which the instruments
rest. \citet{Crawfordetal1998} and \citet{CrawfordWebb2000} have
derived a specific pre-processing scheme based on the recordings of the
pressure component of OBS records and comparison between horizontal and
vertical velocity components. \citet{Nouibatetal2022a} and
\citet{Nouibatetal2022b} applied this pre-processing to the AlpArray
OBS recordings, along with corrections for instrument noise (glitches)
at a few stations. 

These noise reduction procedures were insufficient to ensure the
emergence of clean Rayleigh waves in noise correlations for OBS station
pairs, due to seafloor currents, boat traffic, marine animals and
seismic waves in the water column. Since the recordings of onshore
stations are free of such noise, \citet{Nouibatetal2022a,Nouibatetal2022b}
 have synthesized Rayleigh waves between OBS
stations by using onshore stations as virtual sources. This procedure
is named C2, or iterative noise correlation because it recovers the
Rayleigh-wave signal between two OBSs by correlating Rayleigh-wave
signals emerging from correlations between these OBSs and onland
stations. The C2 method relies on the stationary phase theorem. Onland
stations used as virtual sources must be located close to the azimuth
of the OBS pair to optimize constructive interference of the wavefields
radiated by the source and recorded by the two OBS. In 
\citet{Nouibatetal2022a,Nouibatetal2022b}, virtual sources were
selected in azimuths close (${\pm}$20\textdegree) to the azimuth of the
OBS pair. Since virtual sources are mostly distributed to the North and
East of the OBS array, \citet{Nouibatetal2022b} enhanced the coverage
by separate use of the causal and anticausal parts of the first-order
correlations (C1: between OBSs and land stations) to compute the
OBS--OBS second-order correlations (C2). Extensive tests show that
Rayleigh wave signal quality may be higher in OBS--OBS correlations (C1)
than in C2 correlations in the 5--10~s period range with lower water
column noise. For each OBS--OBS pair, \citet{Nouibatetal2022b} selected
the correlation of highest quality after checking the coherence of C1
and C2 correlations. Figure~\ref{fig2} documents how this procedure
improves inter-OBS signal quality, and thus ray coverage inside the
Ligurian basin, enabling high-resolution ambient noise tomography of
the crust in the basin.


\begin{figure*}
\includegraphics{fig02}
\caption{\label{fig2}Inter-OBS noise correlation waveforms (Z comp.)
generated in four frequency bands by: (left) first-order, standard
correlation of pre-processed noise records; (right) second-order
correlation using land stations as virtual sources. Modified from
\citet{Nouibatetal2022b}.}
\end{figure*}

\subsection{Bayesian and transdimensional inference}\label{sec2.2}

A primary objective of (isotropic) seismic tomography of the Alpine
lithosphere is a depth map of the crust-mantle boundary at a resolution
of a few tens of km. More precisely, a geological
\citep{Calcagnoetal2008} and potentially geodynamic modeling of the
mountain belt requires depth maps of the superposed European, Adriatic
and Ligurian Moho surfaces affected by the subduction of the Alpine
Tethys, the following collision of the European and Adriatic margins,
and the opening of the Ligurian backarc basin. Seismic studies prior to
the AlpArray project mapped the European Moho to a maximum depth of 55
km below the belt axis
\citep[e.g.,][]{Spadaetal2013,Stehlyetal2009,Kastleetal2018}. Receiver
functions of the CIFALPS transect in the south-western Alps provided
the first seismic evidence for continental subduction, including
converted waves on the deep European Moho at 75--80~km beneath the
westernmost Po plain, and negative-polarity conversions on an
``inverted'' Moho between the Adriatic mantle wedge on top and the
European crust below, as shown in Figure~\ref{fig3}b
\citep{Zhaoetal2015,Pauletal2022}. 


\begin{figure*}
{\vspace*{-4pt}}
\includegraphics[scale=0.95]{fig03}
{\vspace*{-5pt}}
\caption{\label{fig3}Cross-sections along the CIFALPS profile in the
southwestern Alps (CC$'$, location shown in Figure~\ref{fig1}b).
(a)~Bouguer anomaly from \citet{Zahorecetal2021}; IGA: Ivrea gravity
anomaly. (b)~Common-conversion point (CCP) stack of receiver functions
at stations of the CIFALPS experiment; locations of stations projected
onto the CIFALPS profile are shown as black inverted triangles at the
surface; the solid thick black lines are the European Moho (EurM) and
the Adriatic Moho (AdM) that appear as converted waves of positive
amplitude (in orange-red); the dotted thick black line shows an
inverted Moho (InvM, higher velocity on top than on the bottom), which
is marked by converted waves of negative amplitude (in blue); modified
from \citet{Pauletal2022}. (c)~Shear-wave velocity model from the
ambient noise tomography of \citet{Luetal2018}. (d)~Shear-wave velocity
model from the ANT with transdimensional inversion to $V_{\mathrm{s}}$ by
\citet{Zhaoetal2020}. (e)~Shear-wave velocity model from the ambient
noise tomography of \citet{Nouibatetal2022a}; dashed black line: $V_{\mathrm{s}}=
3.8$ km/s velocity contour; continuous black line: $V_{\mathrm{s}}= 4.3$~km/s
velocity contour (Moho proxy); the thick grey lines are the Moho
boundaries shown as black lines in (b). (f)~Shear-wave velocity model
from the ambient noise wave-equation tomography of
\citet{Nouibatetal2023}.}
\end{figure*}


To take full advantage of the increased station density and to avoid
dependence on arbitrary choices such as an initial model, several
improvements were made to the inversions for $V_{\mathrm{s}}$ with depth.


\subsubsection{1-D depth inversion for $V_{\mathrm{s}}$ with a full exploration of
the model space}\label{sec2.2.1}

A first improvement of the 1-D inversion for $V_{\mathrm{s}}$ was performed by
\citet{Luetal2018}, who used a subset of AlpArray data, as station
installation was still ongoing. Their inversion for group velocity maps
involved a linearized inversion method based on ray theory
\citep{BoschiDziewonski1999}, and an adaptive parameterization scheme
that reduced cell size in areas with high path density. 

In the 1-D inversion stage for $V_{\mathrm{s}}$, they built upon the inversion method
of \citet{Macquetetal2014}, where a full exploration of the model space
was combined with a linear inversion for $V_{\mathrm{s}}$. In this way,
\citet{Luetal2018}'s inversion for $V_{\mathrm{s}}$ included two steps. The first one
used a grid search approach to uniformly sample and calculate group
velocity dispersion curves for a low-dimensional model space based on a
three-layer crust above a mantle half-space. This full model
exploration was feasible because the dispersion curves can be computed
with normal mode summation, and is computationally tractable. It was
thus possible to determine an ensemble of models for which the
dispersion curve matched the observed one. This first inversion step
provided for each geographical grid cell a probabilistic $V_{\mathrm{s}}$ model and
the probability to have a layer boundary at a given depth. However, to
reduce the size of the model space to explore, the parameterization did
not allow for low velocity zones. The second step was a linear
inversion for $V_{\mathrm{s}}$, where the problem was linearized around the mean of
the probability distribution obtained at the previous step. 

The dense station coverage and the two-step 1-D inversion scheme for $V_{\mathrm{s}}$
combining probabilistic sampling and linear inversion resulted in a
(pseudo) 3-D $V_{\mathrm{s}}$ model with significantly higher spatial resolution
than, for example, the model by \citet{Kastleetal2018}. Comparison with
controlled-source reflection profiles {[}ECORS-CROP:
\citealp{SenechalThouvenot1991}; TRANSALP: \citealp{Kummerowetal2004}]
and receiver function stacked sections across the Alpine belt
{[}TRANSALP: \citealp{Kummerowetal2004}; CIFALPS:
\citealp{Zhaoetal2015,Pauletal2022}{]} displayed striking coincidence
in Moho depth despite the poor sensitivity of surface-wave dispersion
to layer boundaries (see Figures~\ref{fig3}b,c for the CIFALPS
transect). \citet{Luetal2018} computed three depth maps of Moho proxies
in the Alpine region using the probability of interface occurrence, the
depth gradient of $V_{\mathrm{s}}$ and the isovelocity surface $V_{\mathrm{s}}= 4.2$ km/s. These
maps revealed new features such as an 8-km abrupt Moho step beneath the
external crystalline massifs of the Western Alps, from Pelvoux to
Mont-Blanc, which had not been detected by the ECORS-CROP DSS profile
presumably due to poor signal penetration.

The 3-D $V_{\mathrm{s}}$ model by \citet{Luetal2018} however failed to clearly image
the continental subduction of Europe beneath Adria in the western Alps,
which had been identified by receiver functions
(Figures~\ref{fig3}b,c). Since the inversion of \citet{Luetal2018} was
based on a simplified parameterization (only three crustal layers), and
on a final optimization framework, the solution was a unique model that
explains observations, but without associated uncertainties and
trade-offs. In this way, this unique solution did not depict the full
state of information contained in the data. 


\subsubsection{1-D depth inversion for $V_{\mathrm{s}}$ with
transdimensional inference}

The next step in methodological improvement was therefore to carry out
inversions within a full Bayesian framework and a fully adaptive
parameterization, without any linearization, and where the \mbox{solution} is
a full probability distribution, therefore providing uncertainty
estimates. In this context, \citet{Zhaoetal2020} performed a Bayesian
transdimensional inversion of the group-velocity Rayleigh-wave
dispersion data of \citet{Luetal2018} with a focus on the western Alps
(4.5\textdegree{}E--9\textdegree{}E; 44\textdegree{}N--46.7\textdegree{}N). 
Unlike the inversion by \citet{Luetal2018}, which assumed a four-layer
model, the transdimensional inversion treats the number of layers as an
unknown parameter \citep{Bodinetal2012,YuanBodin2018}. At each
gridpoint, \citet{Zhaoetal2020} inverted dispersion data for $V_{\mathrm{s}}$
perturbations around a homogeneous half-space reference model with
velocity 3.8~km/s, allowing a wide range of velocity variations
(${\pm}$50\%). 

The resulting (pseudo) 3-D $V_{\mathrm{s}}$ model (Figure~\ref{fig3}d)
displays a channel of anomalously low shear-wave velocities
($V_{\mathrm{s}}= 3.6$ km/s) at 70-km depth beneath the fast-velocity,
high-density Ivrea body anomaly in the CIFALPS cross-section (IGA in
Figure~\ref{fig3}a). This low-velocity anomaly is located between the
deep European Moho and the inverted Moho of the Ivrea body imaged by
receiver functions of the CIFALPS experiment (Figure~\ref{fig3}b). The
transdimensional inversion of \citet{Zhaoetal2020} was thus able to
image the subduction of the European continental lithosphere beneath
the Adriatic lithosphere with minimal a priori constraints.

\subsubsection{3-D tomography with transdimensional inversions}

The third methodological improvement carried out by
\citet{Nouibatetal2022a} was to use the entire AASN dataset and apply a
full transdimensional approach at both stages of the inversion, that is
in the 2-D inversions for group velocity maps and in the 1-D inversions
for $V_{\mathrm{s}}$. At each period, a transdimensional inversion was carried out to
obtain probabilistic Rayleigh-wave group-velocity maps following
\citet{Bodinetal2012}. In the second stage, a full transdimensional
approach was used to invert for $V_{\mathrm{s}}$ at each geographical location the
dispersion curve obtained in the first stage.

The key benefit of the transdimensional inversion for group velocity
maps is that spatial parameterization is treated as part of the
inversion, allowing local resolution to adapt to data density. Another
improvement introduced by \citet{Nouibatetal2022a} is the computation
of uncertainties on group velocity estimates by transdimensional
inversion; these \mbox{uncertainties} are then incorporated into 
the inversion for $V_{\mathrm{s}}$ (second stage). 
\citet{Youngetal2013} and \citet{Piliaetal2015}
were the first studies to use this two-step transdimensional procedure,
with uncertainty computed in the first stage used to weight the input
in the second.

In the first stage, to account for the strong lateral velocity
contrasts of the Alpine crust, straight rays of the classical forward
model were replaced by bent rays with the ray geometry updated at each
iteration using the fast marching method of
\citet{RawlinsonSambridge2004}. \citet{Nouibatetal2022a} applied this
new inversion scheme to four years of noise records from ${\sim}$1440
permanent and temporary seismic stations, including the entire AlpArray
network with its offshore component (Z3 network: AlpArray Seismic
Network 2015), the two CIFALPS experiments {[}YP and XT networks:
\citealp{Zhaoetal2016b,Zhaoetal2018}], and the EASI experiment {[}XT
network: \citealp{AlpArraySeismicNetwork2014}{]} (Figure~\ref{fig1}).
Data coverage is improved compared to \citet{Luetal2018}, and so to
\citet{Zhaoetal2020}, especially in the Western Alps and the Ligurian
Basin. 

Figure~\ref{fig3}e shows a depth section along the CIFALPS transect in
the 3-D $V_{\mathrm{s}}$ model by \citet{Nouibatetal2022a},
which may be compared to
the models by \citet{Luetal2018} shown in Figure~\ref{fig3}c, and 
\citet{Zhaoetal2020} shown in Figure~\ref{fig3}d. Unlike
\citet{Luetal2018} and similar to \citet{Zhaoetal2020,Nouibatetal2022a}
image the dipping low-velocity layer beneath the Ivrea body
high-density, high-velocity anomaly (IGA, $x=210$--300~km in
Figure~\ref{fig3}), which is indicative of the continental subduction
of Europe beneath Adria. As the \citet{Zhaoetal2020} inversion for $V_{\mathrm{s}}$
was based on the same dispersion dataset as \citet{Luetal2018}, and as
\citet{Nouibatetal2022a} used the same inversion method for $V_{\mathrm{s}}$ as
\citet{Luetal2018}, but with a different set of four-layer models, we
propose that the major difference between Figures~\ref{fig3}c
and~\ref{fig3}e is related to the difference between the sets of
velocity models explored in the probabilistic inversion for $V_{\mathrm{s}}$: 
130~million for \citet{Nouibatetal2022a}, and 8~million for
\citet{Luetal2018}. In particular, the complex vertical velocity
profiles of the subduction region with alternating high and low
velocities were not explored by \citet{Luetal2018}. Vertical velocity
gradients are stronger in \citet{Nouibatetal2022a}'s model
(Figure~\ref{fig3}e) than in \citet{Zhaoetal2020}'s model
(Figure~\ref{fig3}d), in particular at the crust-mantle boundary. Even
though the input data and the inversion schemes differ, the overall
similarities between the two models of \mbox{Figures~\ref{fig3}e 
and~\ref{fig3}f} suggest that
the differences in velocity gradient at the Moho may result from the
parameterization of the inversions for $V_{\mathrm{s}}$. While the probabilistic
inversion of \citet{Nouibatetal2022a} assumes a four-layer starting
model with a velocity jump at Moho, the reference velocity model in the
transdimensional inversion of \citet{Zhaoetal2020} is a homogeneous
half-space with no discontinuity. This a priori uniform distribution of
velocity may favor a smooth velocity gradient at the Moho. The $V_{\mathrm{s}}$ model by
\citet{Nouibatetal2022a} with its a priori imposed sharp
crust-to-mantle velocity contrast is more in line with previous
controlled-source (ECORS-CROP in the north-western Alps) and receiver
function profiles [CIFALPS: Figure~\ref{fig3}b; CIFALPS-2,
\citealp{Pauletal2022}], which display clear reflected and converted
phases on the Moho. The impact of other refinements in the
\citet{Nouibatetal2022a} inversion methodology (error estimates on
group velocities, bent rays) is difficult to evaluate due to the
differences in the input data coverage. It would be an interesting test
to apply the fully transdimensional inversion for $V_{\mathrm{s}}$ of
\citet{Zhaoetal2020}, with an a priori constraint on the velocity jump
at the Moho, to the group velocity data with uncertainties of
\citet{Nouibatetal2022a}.


\subsection{Wave-equation tomography}

To compute shear-wave velocity models of the crust and uppermost mantle
from Rayleigh-wave traveltime data between station pairs,
\citet{Luetal2018} and \citet{Nouibatetal2022a} have used a two-stage
inversion scheme, which, in spite of improvements related to a
probabilistic approach, is a standard strategy for ANT. The first
stage, a series of 2-D inversions of traveltime data for group velocity
maps at selected periods, is based on ray theory, which is only valid
at infinite frequency. The second stage, a set of 1-D inversions of the
local dispersion curves for S-wave velocity, results in 1-D models
merged into a pseudo 3-D model. The strong crustal heterogeneity of the
Alps and surrounding areas, including the Ligurian back-arc basin,
makes the ray hypothesis particularly inadequate. This area therefore
provides an excellent chance of testing a tomography procedure that\break 
better accounts for the physics of seismic wave\break propagation. In this
section, we first describe a new wave-equation based approach of
ambient noise tomography called wave equation tomography {[}WET,
\citealp{Luetal2020}{]} for the elastic case, using only land station
recordings. We next present the extension of the WET method to include
elastic--acoustic coupling at the sea bottom and OBS records
\citep{Nouibatetal2023}.

\subsubsection{Wave-equation onshore tomography (elastic case)}

To calculate a truly three-dimensional $V_{\mathrm{s}}$ model consistent with wave
propagation physics, \citet{Luetal2020} derived a wave-equation based
approach, hence called ambient noise wave-equation tomography. As in
the ambient noise adjoint tomography of \citet{Chenetal2014}, the
observable is travel time (i.e.\ phase) of the Rayleigh wave
reconstructed from noise correlation, while full-waveform inversion
(FWI) would invert amplitude and time. Indeed, Rayleigh wave amplitudes
are not correctly retrieved by the classical noise correlation
procedure \citep[e.g.,][]{Campillo2006}. FWI of ambient noise correlation
signals probably has great potential for lithospheric imaging. It is
not yet operational because it would require specific pre-processing of
noise records for a better retrieval of amplitudes and, more
importantly, an accurate estimate of noise source distributions and
emitted waveforms \citep[e.g.,][]{Fichtner2014,Sageretal2018}. 

The wave-equation tomography (WET) approach implemented by
\citet{Luetal2020} consists of iteratively updating an initial ambient
noise tomography model (from the two-stage traditional method) by
minimising frequency-dependent phase traveltime differences between the
Rayleigh waveforms observed in noise correlations and the synthetic
waveforms computed by 3-D numerical modeling of wave propagation.
Synthetic waveforms are computed with the 3-D elastic wave equation
solver of the SEM46 package \citep{Trinhetal2019} based on the spectral
element method \citep{KomatitschVilotte1998,KomatitschTromp1999}. Of
the 304 stations located in the study region, a subset of 64 suitably
located stations were selected as virtual sources by
\citet{Luetal2020}, and the signal bandwidth was limited to
{[}10--50~s{]} to ensure acceptable computing time. The observed signal
for a station pair is the Rayleigh wave part of the Green's function,
estimated from the time derivative of the cross-correlation of seismic
noise. The synthetic signal is the convolution product of a bandpass
filtered Dirac delta source function with the \mbox{synthetic} Green's
function computed with SEM46 for a source--receiver pair, and a vertical
force applied on the free surface at the source location. The misfit
function is then computed from the frequency-dependent phase traveltime
differences between observed and synthetic waveforms. The adjoint-state
approach is used to compute the misfit gradients
\citep[e.g.,][]{Trompetal2005}, and the inversion is conducted as an
iterative local optimization problem.

The final S-wave velocity model obtained by \citet{Luetal2020} after 15
iterations of wave equation tomography has a 65\% lower total misfit
than the initial ANT model by \citet{Luetal2018}. This strong misfit
reduction is mostly due to periods larger than 25~s, where the WET has
corrected for a strong and unexplained positive shift of the traveltime
misfit histograms in the initial model. A direct consequence is that
the final $V_{\mathrm{s}}$ model has significantly higher average velocities at lower
crustal depth (30 km). Beyond this tuning of average velocities, the
WET approach was able to retrieve finer scale heterogeneities than the
ANT model, for example a high-velocity anomaly at 10 km depth, which is
closer in shape to the well-known Ivrea positive Bouguer anomaly 
\citep[see Figure~8 in][]{Luetal2020}.

\subsubsection{Wave-equation onshore/offshore tomography (with
acoustic--elastic coupling)}

As compared to \citet{Luetal2018}, a more accurate ANT model of the
Alpine region and its surroundings was computed by 
\citet{Nouibatetal2022a} using a larger noise correlation dataset including
recordings of sea-bottom seismometers in the \mbox{Ligurian} basin, and an
improved inversion scheme (see Section~\ref{sec2.2}). Like
\citet{Luetal2020}, \citet{Nouibatetal2023} have performed a WET to
improve the model of \citet{Nouibatetal2022a} by accounting for
three-dimensional and finite-frequency effects of wave propagation. A
major improvement achieved by \citet{Nouibatetal2023} is the inclusion
of the water layer effect on wave propagation in the Ligurian basin.
The vertical-component record of short-period (${<}$20~s) surface waves
at ocean-bottom stations is dominated by a fluid--solid interface wave
named Rayleigh--Scholte wave that propagates at lower speed than the
Rayleigh wave in the elastic medium. At periods ${>}$20~s, the effect of
the water layer becomes negligible for the Ligurian basin (max.\ water
depth 3 km) because the Rayleigh wavelength is large compared to water
depth \citep{Nouibatetal2023}. Accounting for the effect of the water
layer in the forward simulation of surface-wave propagation and in the
inversion of traveltime observations for shear-wave velocity is
therefore required to enhance resolution in the shallow layers of the
Ligurian basin.

For their wave-equation tomography with elastic--acoustic coupling,
\citet{Nouibatetal2023} selected, from the dataset of
\citet{Nouibatetal2022a}, 600 stations as receivers, out of which 185
were selected as sources based on noise correlation signal quality. The
initial model was the ANT $V_{\mathrm{s}}$ model computed by
\citet{Nouibatetal2022a}. As in \citet{Luetal2020}, the inverse problem
minimised the frequency-dependent phase traveltime differences between
observed and synthetic vertical-component waveforms for considered
source--receiver pairs. The inversion was conducted in the 5--85~s period
range, considering elastic--acoustic coupling in the forward simulation
in the 5--20~s band. Unlike \citet{Luetal2020}, the velocity model was
updated progressively from long periods (40--85~s) to shorter ones,
20--40~s, 10--20~s then 5--10~s to avoid cycle skipping. 

The final WET $V_{\mathrm{s}}$ model differs significantly from the initial ANT $V_{\mathrm{s}}$
model, particularly in shallow layers of the Ligurian Basin, and in the
most heterogeneous parts of the Alpine crust. These discrepancies
highlight the importance of accounting for the physics of wave
propagation, i.e.\ elastic--acoustic coupling at the sea-bottom,
finite-frequency effects and 3-D propagation. In the alpine crust (i.e.\ 
onland), velocity contrasts tend to be slightly enhanced, in particular
at lower crustal depth [26 km, see\break Figure~6 in
\citet{Nouibatetal2023}], while the location and shape of velocity
changes are preserved with respect to the initial ANT model
(Figures~\ref{fig3}e--f). At depths shallower than ${\sim}$10 km, S-wave
velocities in the Ligurian basin are ${\sim}$8\% lower in the WET model
than in the ANT model, in agreement with the lower velocities of the
Rayleigh--Scholte wave. Accounting for the water layer in the Ligurian
basin also leads to higher velocities in the crust of western
(Variscan) Corsica.

\citet{Nouibatetal2023} conducted extensive evaluation of the
robustness of their WET model. The quality of the S-wave velocity model
was documented by traveltime misfit maps at representative periods
between 8 and 55~s. As expected, misfit is lower for the WET model than
for the initial ANT model, except in the peripheral poorly illuminated
regions \citep[see Figure~12 in][]{Nouibatetal2023}. The footprint of
the broad geological structure has disappeared, whereas it was clearly
visible for the ANT model at periods ${\geq}$20~s. Using the weak
sensitivity of Rayleigh wave phase velocity to $V_{\mathrm{p}}$,
\citet{Nouibatetal2023} also inverted Rayleigh wave dispersion
observations for P-wave velocity, starting from an initial $V_{\mathrm{p}}$ model
computed from the initial (ANT) $V_{\mathrm{s}}$ model using an empirical formula. To
document the robustness of their P and S-wave models, they computed
synthetic waveforms for a regional earthquake in Switzerland and
compared to observed waveforms in three frequency bands between 2 and
10~s. The fit between simulated and observed seismograms is striking for
the travel times of the P and Rayleigh waves, but also for relative
amplitudes between the wave trains. Such a high coherence in the
short-period bands is a reliable indication that wave equation
tomography is capable of imaging small-scale crustal structures. This
synthetic test also demonstrates that the WET $V_{\mathrm{s}}$ and $V_{\mathrm{p}}$ models derived
from ambient noise correlations would be a good initial model for a
wave-equation tomography using earthquake records. As regional
earthquake records include both body and surface waves, the resulting
P-wave velocity model would be much more accurate than that obtained
from Rayleigh wave travel time inversion alone, due to its weak
sensitivity to $V_{\mathrm{p}}$.

\subsection{Model validation: Moho depth}\label{sec2.4}

Surface waves are weakly sensitive to velocity contrasts, but Bayesian
inversions for $V_{\mathrm{s}}$ were shown to be very effective for imaging the
crust-mantle boundary \citep{Luetal2018,Nouibatetal2022a}.
Comparisons of Moho depths imaged by other geophysical methods such as
DSS profiles (ECORS-CROP) and receiver function analyses have shown
that selected iso-velocity surfaces in the ambient noise derived 3-D $V_{\mathrm{s}}$
models are reliable proxies of the Moho
\citep{Nouibatetal2022a,Pauletal2022,Nouibatetal2023}. For example, the
velocity contour $V_{\mathrm{s}}= 4.3$~km/s is a good proxy for the European Moho,
and the Adriatic Moho outside the Ivrea body\break region, as shown by the
coincidence of this surface with Moho conversions picked from CCP
receiver\break function stacks along the CIFALPS and CIFALPS2 sections (see
for example Figure~\ref{fig3}e). In the Ivrea body region, the Moho
imaged by the ECORS-CROP DSS profile and receiver function sections is
better approximated by the 3.8~km/s velocity surface, probably because
the external rim of the Ivrea body is made of serpentinized peridotites
with lower velocities than dry peridotites {[}Figure~\ref{fig3}e; see
discussion in \citealp{Malusaetal2021}{]}. In the Ligurian basin, a
comparison of the $V_{\mathrm{s}}$ model of \citet{Nouibatetal2022a} with a $V_{\mathrm{p}}$ model
derived by \citet{Dannowskietal2020} from refraction and wide-angle
reflection OBS data shows a striking coincidence of the
$V_{\mathrm{p}}$ Moho (iso-velocity 7.2~km/s) with the 4.1 km/s
$V_{\mathrm{s}}$ contour  \citep[see Figure~9 in][]{Nouibatetal2022b}.
A Moho depth map of much higher resolution than previous ones
\citep[e.g.,][]{Gradetal2009,Spadaetal2013} can therefore be built by
mapping iso-velocity contours of the $V_{\mathrm{s}}$ models by
\citet{Nouibatetal2022a} and \citet{Nouibatetal2023}. 

\section{Anisotropic ambient noise tomography}\label{sec3}

Surface waves in seismic noise correlations can potentially shed new
light on anisotropy in the Earth's crust and upper mantle, overcoming
the observational biases related to the uneven distribution of
earthquakes, which leads to uneven azimuthal coverage 
\citep[for a review, see][]{MaupinPark2015}. 

Our knowledge of azimuthal anisotropy beneath the Alps and Italy mostly
comes from XKS splitting studies. Fast velocity directions are
generally parallel to the mountain chain in the Western and Central
Alps
\citep{Barruoletal2004,Barruoletal2011,Lucenteetal2006,Salimbenietal2018}, 
the Eastern Alps \citep{Bokelmannetal2013,Qorbanietal2015} and the
Apennines \citep{Palano2015}. Teleseismic P wave travel time delays
were used by \citet{Rappisietal2022} to study anisotropy in the Central
Mediterranean area, including the Alps. Both XKS and teleseismic P
waves have too steep incidence angles to reliably constrain anisotropy
in the crust, also because the cumulated travel time in the crust is
much smaller than in the upper mantle. Studies using regional refracted
Pn and Sn phases like \citet{Diazetal2013} provide information about
the uppermost mantle, as these waves are tied to the Moho interface and
propagate with mantle velocities. 

Long wavelength anisotropic structure beneath Europe is mostly known
from large scale or global surface wave studies where the crust is not
resolved
\citep[e.g.,][]{Kustowskietal2008,WeidleMaupin2008,Boschietal2009,Zhuetal2015,Nitaetal2016,SchivardiMorelli2011}. 
All these studies are constructed from long period observations, where
the crust is not inverted for but set from a reference model such as
LITHO1.0 \citep{Pasyanosetal2014}. A few studies focussing on crustal
anisotropy used surface waves from seismic noise correlations to
characterise anisotropy at regional scale, such as Switzerland
\citep{Fryetal2010}, the Vienna Basin \citep{Schippkusetal2020}, the
Bohemian Massif \citep{Kvapiletal2021}, or the Eastern Alps
\citep{Kastleetal2024}. 



The dense station coverage (Figure~\ref{fig1}) across the greater
Alpine region during the AlpArray project provided a unique opportunity
to further develop and understand the limitations of noise-based
imaging methods aimed at characterising seismic anisotropy across a
strongly heterogeneous structure. Since the Alpine region is
heterogeneous at all scales, it is in particular crucial to reliably
estimate uncertainties and identify locations and periods where
systematic errors may affect azimuthal anisotropy measurements. These
questions are the focus of Section~\ref{sec3.1}. 

Once uncertainties associated with surface wave velocities are
estimated, they can be used in inversions for depth variations of
anisotropic parameters. Reliable uncertainties are also crucial in the
case of joint inversions, where different datasets are inverted
simultaneously (e.g., earthquake based and noise-based observations). A
Bayesian inversion framework in theory handles this question in a both
intuitive and simple way, producing families of acceptable Earth
models, over which it is possible to make posterior statistics. This
ensemble solution can be exploited to give reliable and useful
information on the Earth structure. Section~\ref{sec3.2} is dedicated
to such strategies for both radial and azimuthal anisotropy.

\subsection{Improving observations of azimuthal\newline
anisotropy}\label{sec3.1}

The effect of isotropic heterogeneities on the wavefield can map into
anisotropic parameters. In this section, we give some insights into
such effects through the now well-established Eikonal tomography method
\citep{Linetal2009}. We present outcomes of Eikonal tomography across
the Alps and show how beamforming can improve seismic noise-based
observations of anisotropy. Both Eikonal \mbox{tomography} and beamforming
naturally correct for deviations from great circle
\citep[e.g.,][]{Pedersenetal2015} and uneven distribution of noise sources
\citep[e.g.,][]{Fromentetal2010,Harmonetal2010}, so the only main bias to
handle is the one from isotropic heterogeneities.

\vspace{0.5em}
\subsubsection{Biases in measurements of azimuthal\newline anisotropy:
insights from Eikonal tomography}\label{sec3.1.1}

Eikonal tomography makes use of dense station networks to reconstruct
the wavefield propagating away from a source \citep{Linetal2009}. In
ambient noise applications, any seismic station can act as a\break virtual
source and the signals recorded at all other stations allow us to image
the wavefield. When computing ambient noise correlations, amplitudes
are lost due to filters such as spectral whitening. Only the
travel-time field is recovered, therefore hindering\break the use of
amplitudes to correct for bias from wave interference. 

In Eikonal tomography \citep{Linetal2009}, travel times recorded at
receivers are interpolated onto a regular grid. The gradient of the
travel-time field then gives the phase velocity and the propagation
direction at each grid point. This highlights important strengths of
Eikonal tomography: it is simple to calculate and the direction of the
incoming wave is directly determined from the data. In order to get the
full, azimuthal anisotropic phase velocity map, the process has to be
repeated for all available virtual sources so that a set of phase
velocity measurements and their propagation azimuths is obtained for
each grid cell. The final phase-velocity field is obtained by averaging
the measured phase velocities for all virtual sources. Different
statistics of the estimated phase velocities (mean, median, mode and
standard deviation) can be used as uncertainty estimates. Similarly,
for the anisotropic part, the residual from the fitted function
provides information on the reliability of parameter estimates.

Isotropic bias can be approximately modelled by a 360\textdegree\ 
variation of phase velocity with azimuth
\citep{LinRitzwoller2011,Mauerbergeretal2021}. We refer to phase
velocity variations with azimuth in the following form (applicable to
weakly anisotropic media):
{\begin{eqnarray}
C(\phi )&=&C_{0}\left(1+A_{1}\cos(\phi -\theta _{1})+A_{2}\cos(2(\phi -\theta
_{2}))\right.\nonumber\\
&&\left.+\,A_{4}\cos(4(\phi -\theta _{4}))\right),\label{eq1}
\end{eqnarray}}\unskip
where $C(\phi )$ {is the anisotropic phase velocity, dependent on the
wave propagation azimuth} $\phi $. It is given by its isotropic
component $C_{0}$ {plus three anisotropy terms that describe the}
$\theta _{1}$ {anisotropy (360\textdegree~symmetry), the} $\theta _{2}$
{anisotropy (180\textdegree~symmetry) and the} $\theta _{4}$ {anisotropy
(90\textdegree~symmetry), where} $\theta $ {is the fast axis direction
and} $A$ {the fast axis amplitude in fractions of} $C_{0}$ {from zero
to peak. The} $\theta _{2}$ {anisotropy is the dominant component for
Rayleigh anisotropy \citep[e.g.,][]{MontagnerNataf1986} while the
non-physical} $\theta _{1}$ anisotropy can be used to check for
potential measurement bias created by isotropic heterogeneities. The
$\theta _{4}$ {component is not included in further discussions since
it is mostly relevant for Love waves which are not taken into account
in the methods described below.}\looseness=-1


Figure~\ref{fig4} shows systematic errors resulting from isotropic
heterogeneities through synthetic data computed in a laterally
heterogeneous model, onto which we apply Eikonal tomography. In
Figure~\ref{fig4}d (perfect sampling), we observe an undulating pattern
in the isotropic velocities, which depends on the signal wavelength.
This undulating pattern is due to the interference of waves that are
refracted by the isotropic anomaly. Most of the bias, however, cancels
out when averaging over sources from many directions
(Figure~\ref{fig4}e). This error can also be corrected when taking the
amplitude information into account 
\citep[e.g., Helmholtz tomography,][]{LinRitzwoller2011}.
A more severe bias stems from
under-sampling the wavefield (Figure~\ref{fig4}b). To work properly, the
Eikonal method requires approximately plane wavefront between adjacent
receivers. If the medium is highly heterogeneous, this assumption is
violated and only a smoothed version of the true wavefield can be
reconstructed from the measurements at the receiver locations. The
resulting error in the isotropic part mostly cancels out when signals
from different azimuths are averaged (Figure~\ref{fig4}c). In the
anisotropic part, however, a systematic, azimuth-dependent error
remains, resulting in spurious anisotropy around the edges of strong
velocity heterogeneities.
%\looseness=-1


\begin{figure*}
\vspace*{-5pt}
\includegraphics{fig04}
\caption{\label{fig4}Example illustrating the bias in the anisotropic
Eikonal tomography method. (a)~2-D wavefield propagating away from the
virtual source station; black dashed lines indicate the locations of
isotropic high and low velocity anomalies. (b)~Isotropic phase
velocities recovered by applying the Eikonal tomography method to the
wavefield recorded at the receiver locations; anomalies outside the
dashed lines are mainly caused by incomplete reconstruction of the true
traveltime field. (d)~Same as (b) for a perfectly sampled traveltime
field; the striped velocity heterogeneities behind the anomalies are
caused by wavelength-dependent wave interference effects. (c) and (e)
Final models resulting from averaging the models for all virtual
sources, i.e., swapping the position of the virtual source with all
available receiver positions. The apparent anisotropy (yellow bars) is
caused by velocity errors illustrated in (b) and (d) since the input
model is purely isotropic.}
\end{figure*}


The level of these spurious anisotropic measurements depends not only
on the absolute velocity variation but also on the wavelength of the
signal, the inter-station spacing and the velocity gradient at the
edges of the anomalies. The controlling factor for the bias is the
level of complexity of the wavefront between adjacent stations.

\subsubsection{Eikonal tomography in the Alps}\label{sec3.1.2}

\citet{Kastleetal2022} presented an application of Eikonal tomography
to the Alpine area, including data from the AlpArray network. They were
able to resolve both the isotropic and the anisotropic parts at periods
sensitive to different depth ranges from the shallow crust to the
uppermost mantle. The three anisotropic components were obtained by
taking the azimuthally dependent phase-velocity measurements and
fitting them to Equation~(\ref{eq1}). 


\citet{Kastleetal2022} tested for different observables that can be
related to a bias in anisotropy measurements, such as the amplitude of
the $\theta _{1}$ component or the standard deviation from fitting
the $\theta _{2}$ {anisotropy. They found that none of these
observations are adequate to reliably identify biased measurements.
More sophisticated approaches such as forward modelling the travel-time
field from the isotropic velocities may be necessary. This could also
be used iteratively to improve the accuracy of the interpolation for
highly complex travel-time fields, at the cost of giving up the
relative simplicity of the Eikonal approach. However,
\citet{Kastleetal2022} inferred from synthetic tests that isotropic
anomalies of up to 10\% produce bias, which in most cases is smaller than
what is identified as ``true'' anisotropy.} 


{The} $\theta _{2}$ anisotropic component obtained by
\citet{Kastleetal2022} is presented in Figure~\ref{fig5}. The fast axis
directions generally follow the curvature of the Alpine arc at 15~s
period, representative of the upper crust. At 40~s period, the fast
axis becomes more oblique to orogen-perpendicular within the central
Alps. The eastern part of the Alps has more homogeneously E--W oriented
fast axes from short to long period measurements. In most of the
northern Alpine foreland, the anisotropic fast axis orientations are
consistently following the Alpine curvature at periods between 10 and
50~s.


\begin{figure*}
\includegraphics{fig05}
\caption{\label{fig5}Maps of $\theta _{2}$ {azimuthal anisotropy at 15
and 40 s period resulting from Eikonal tomography in the Alps. The
direction of each line shows the fast direction, and the length shows
the amplitude (zero to peak). The colour of the lines gives the
estimated uncertainty on the amplitude} $A_{2}$ {(in percent of}
$C_{0}$).}
\end{figure*}


The uncertainties in the anisotropic measurements are largest at model
boundaries, where the azimuthal coverage is non-optimal, as well as
within and around the large sedimentary Po and Molasse basins north and
south of the Alps. In the former, isotropic velocities indicate that
velocities can be as low as ${-}$30\% compared to the average velocity
\citep{Kastleetal2022}. These regions are therefore prone to
the systematic bias mentioned above, particularly at short periods. In
most of the other regions and at longer periods, the anomaly strength
is usually within ${\pm}$10\% with smoothly varying isotropic velocities.
Correspondingly, \citet{Kastleetal2022} found that the imaged
pattern of anisotropies in the Alps is compatible with previous, more
localised studies \citet{Fryetal2010},
but also with the results from beamforming [\citealp{Soergeletal2023}; 
Section~\ref{sec3.1.3}{]}.

\subsubsection{Beamforming}\label{sec3.1.3}

An alternative approach to Eikonal tomography is to use the phase
information in a beamforming approach. Beamforming consists of, for a
set of closely located seismic stations and at a given frequency,
inferring the direction and local phase velocity of the incoming wave
\citep[for a review, see][]{RostThomas2002}]. Contrary to Eikonal
tomography, the full waveform of the empirical Green function is used,
as illustrated in Figure~\ref{fig6}. The original waveforms\break
(Figure~\ref{fig6}a) are shifted in time using a grid search across
many phase velocities and incoming directions. The final estimated
phase velocity is 
\begin{figure*}%[H]
{\vspace*{5pt}}
\includegraphics{fig06}
{\vspace*{5pt}}
\caption{\label{fig6}Illustration of data processing, and map of
estimated azimuthal anisotropy at 15~s period for one seismic array.
(a)~Filtered traces before aligning. (b)~Same traces as (a), aligned
using the combination of azimuth and phase velocity for which the stack
is optimal (maximum peak amplitude). (c)~Stack at optimal parameters,
normalised by the number of stations. The maximum of this stack is
called $d_{\max}$. (d)~Phase velocity as a function of azimuth, using
all source stations. The orange and green curves correspond to the
first two cosine terms of Equation~(\ref{eq1}). (e)~Number of source
stations (binned in 5\textdegree~azimuth bins) as a function of azimuth.
(f)~Quality of stack $d_{\max}$ for all source stations. The seismic
array is centered at the blue dot. It is surrounded by a blank area,
because of the need to respect a minimum distance between source and
receiver stations. The colour code shows $d_{\max}$, the maximum
amplitude of the stack divided by the number of traces.
\citet{Soergeletal2023} particularly focused on quality control and
error estimates to be able to carry out inversions for anisotropic
parameters with depth. Such quality control and error estimates are key
input for subsequent transdimensional Bayesian inversion described in
Section~\ref{sec3.2}.}
{\vspace*{6pt}}
\end{figure*}
the one for which the traces are optimally aligned
(Figure~\ref{fig6}b), and for which the stack yields the highest
maximum amplitude (Figure~\ref{fig6}c). 


Beamforming has been extensively used on seismic noise, in particular
to characterise the noise field and infer the locations of origin of
noise sources {[}for Europe, see for example
\citealp{Friedrichetal1998,Landesetal2010} or
\citealp{JuretzekHadziioannou2016}{]}. Beamforming has also been used
to quantify anisotropy on earthquake data
\citep[e.g.,][]{Pedersenetal2006,AlvizuriTanimoto2011} or directly on raw
noise records \citep[e.g.,][]{RiahiSaenger2014}. 

To estimate surface wave azimuthal anisotropy, beamforming has to be
carried out, for each source, across the range of target periods
(Figure~\ref{fig6}d).


This method was previously used on earthquake data. Applying it to
noise correlation is straightforward as seismic stations can be used as
virtual sources. Here we present key elements of the implementation of
\citet{Soergeletal2023} and the application to the Alps. A somewhat
similar approach was taken at the same time by \citet{Wuetal2023} for
the Northeastern Tibetan Plateau. 



The specific choices made to adapt the beamforming included: 
\begin{itemize}
\item Using seismic stations outside the study area as virtual sources,
to improve azimuthal coverage at the edge of the study area.

\item Taking into account the short distances between the array and the
``source'' station. This was done mainly by assuming circular wavefronts
following \citet{Maupin2011} rather than plane wavefronts, and imposing
a minimum distance between array and ``source'' station. The
implementation still allowed to take into account great-circle
deviations, which can be significant across AlpArray.

\item Propagating uncertainty of each data point (phase velocity at a
given period for a given array and source stations, see example in
Figure~\ref{fig6}f). This is done by weighting each observed phase
velocity for the estimate of the $A_{2}$ and $\theta _{2}$ terms for a
given array at a given period, with the normalised beam amplitude
$d_{\max}$, and use of bootstrap to estimate the uncertainty on $A_{2}$
and $\theta _{2}$.

\item Exclusion of data points if the isotropic bias (estimated by
$A_{1}$) is significant.

\item Adapting array geometry to wavelength and noise levels.
\end{itemize}


Figure~\ref{fig7} shows two examples of surface wave azimuthal
anisotropy measured with this method across the Alpine area, at 15 s
and 40 s period. Points with risk of a strong isotropic bias (high
value of $A_{1}$ as compared to $A_{2}${, see Equation~(\ref{eq1})) are
excluded from the maps. In many points, the bias was particularly
strong at periods which are sensitive to changes in Moho depth, and was
overall coherent with areas of strong isotropic velocity gradients.}
For further discussion, we refer to \citet{Soergeletal2023}. 



\begin{figure*}
\includegraphics{fig07}
\caption{\label{fig7}Maps of azimuthal anisotropy at 15 s and 40 s
period resulting from the beamforming method. The direction of each
line shows the fast direction, and the length shows the amplitude. The
colour of the lines show the estimated uncertainty on the amplitude
$A_{2}$ (in percent of $C_{0}$).}
\end{figure*}

The beamforming yields spatially coherent patterns of azimuthal
anisotropy across all of the study area. Overall, the beamforming is in
good agreement with the Eikonal method (Figure~\ref{fig5}) on key
findings, such as (at 15~s period) lower amplitude anisotropy within
the Alps as compared to surrounding areas, and the surrounding areas
being dominated by chain parallel fast direction. There are main
differences at the periphery of the study area, which can be explained
by the lack of azimuthal coverage in the Eikonal tomography. Lateral
smoothing effects are also different between the two methods. While the
reliability of the beamforming is high due to the point by point
(geographically, for each period) evaluation of quality and errors, it
also leads to changing geographical distributions of reliable data with
period, and the implementation would need to resolve this issue to be
able to provide 3-D models of anisotropy. 


{Depth inversions for azimuthal anisotropy must be considered a priority
in anisotropic studies because of the wide range of possible values of}
$A_{2}$. For example, if the upper crust is strongly anisotropic and
neither lower crust or upper mantle are anisotropic, the observed
anisotropy at, for example, 40 s period still has its origin in the
upper crust. Similarly, changes in fast direction in layers with small
anisotropy are not directly visible at any given period. Whether small
amplitude anisotropy is resolved therefore strongly depends on the
observation error. The strong control of errors provided by beamforming
is consequently key for the depth inversions which are addressed in the
following section. 
{\vspace*{-7pt}}

\subsection{Transdimensional Bayesian inversions of anisotropic
measurements}\label{sec3.2}

Similar to isotropic inversions, once surface wave velocities have been
estimated at different periods, they can be inverted to recover
anisotropic shear wave velocity structure at depth. This is usually
done in a second inversion step carried out at each geographical
location. Rayleigh and Love dispersion curves can be jointly inverted
for radial anisotropy parameters with depth. In the case of azimuthal
anisotropy, Rayleigh and Love wave dispersion curves and their
azimuthal dependency at each period can each be inverted to obtain
azimuthal anisotropy parameters with depth. 

Fortunately, this inverse problem is only weakly nonlinear, and many
linearized approaches have been used, where data derivatives
(sensitivity kernels) computed around a reference model are used.
However, the solution is highly non-unique, as surface waves are only
sensitive to integrated\break velocities along a depth range. In the case
where seismic anisotropy is inverted for, the non-uniqueness of the
solution becomes particularly problematic as strong trade-offs emerge
between different parameters. For example, the level of radial
anisotropy trades off with $V_{\mathrm{p}}$ and with the level of heterogeneities in
$V_{\mathrm{s}}$ \citep{Bodinetal2015,Alderetal2017,GaoLekic2018}, because of the
equivalence of a stack of horizontal layers and a homogeneous
anisotropic medium \citep{Backus1962}. 

To add extra constraints on the solution, some regularisation can be
used at the cost of biasing estimated uncertainties in the recovered
model. In this way, linearized approaches only provide a unique
velocity model that does not represent the range of potential
solutions, and that is strongly dependent on the regularisation, the
parameterization, or the reference model. For example, the $V_{\mathrm{p}}$ model is
often set to a reference or arbitrarily scaled to $V_{\mathrm{s}}$ variations. The
level of smoothness is also set in advance. Although these choices are
based on geological or mineralogical arguments, they make the
interpretation of results dependent on a priori choices.

Similar to isotropic inversions, these issues can be addressed by fully
exploring the range of potential solutions with Monte Carlo methods
(see Section~\ref{sec2.2.1}). Since forward simulations are
computationally cheap, a large number of 1-D models with variable
parameterizations can be tested to sample complex posterior
distributions. Adaptive parameterizations can be used to explore
complex trade-offs between model parameters. For example, a
transdimensional parameterization (where the number of model parameters
is variable) can be used to explore the ambiguity between the level of
spatial heterogeneity (horizontal layering) and the level of anisotropy
created by the orientation of anisotropic minerals 
\citep{Bodinetal2015,Bodinetal2016,Alderetal2017}. In this case, the
total number of layers in the inverted model as well as the presence of
anisotropy in each layer are not constant parameters. Instead, they are
adjusted by the inversion to fit the data to the degree required by
their estimated noise. In this context, the Bayesian framework enables
to propagate estimated uncertainties in the observed dispersion curves
towards uncertainties in shear wave velocities at depth.


\subsubsection{Adaptive parameterizations: example with radial
anisotropy}\label{sec3.2.1}

To illustrate the benefit of such flexible parameterizations, we first
show in Figure~\ref{fig8} an example of a synthetic test taken from
\citet{Alderetal2021}. A target model is designed with a radially
anisotropic layer in the crust and an isotropic upper crust and mantle.
Noisy synthetic data are created and inverted with a transdimensional
Monte Carlo algorithm, but with two different types of
parameterizations (left and right panels). On the left, the number of
layers is variable and each layer can be either isotropic and described
solely by $V_{\mathrm{s}}$ and $V_{\mathrm{p}}$ (in this case, radial
anisotropy in the form $V_{\mathrm{SH}}/V_{\mathrm{SV}}$ takes a
constant value of 1), or radially anisotropic and described by three
parameters: $V_{\mathrm{SV}}$, $V_{\mathrm{p}}$ and the level of radial
anisotropy $V_{\mathrm{SH}}/V_{\mathrm{SV}}$. In this way, the number
of inverted physical parameters in each layer is variable. On the right
panel, the number of layers is also variable but all layers are
radially anisotropic and described by three parameters. Everything else
is equal in the two inversions.{\looseness=-1}  


\begin{figure*}
\vspace*{-1pt}
\includegraphics{fig08}
\caption{\label{fig8}Joint inversion of Love and Rayleigh wave
synthetic dispersion curves. The same data produced by the true model
are inverted with two different procedures. Posterior distributions of
$V_{\mathrm{SV}}$ (red), and $V_{\mathrm{SH}}/V_{\mathrm{SV}}$
(blue). The true model used to create synthetic data is the thick solid
line in every panel and the mean of the prior distribution used in the
inversion is the grey line in the $V_{\mathrm{SV}}$ panel. Posterior
distributions are depicted with their median (thin solid line) and
likelihood intervals: for each parameter, the dark surface includes
65\% of the models in the ensemble solution while the light area
includes 95\% of the models. Left panels: inversion used in 
\citet{Alderetal2021}, where each layer can either be isotropic or
anisotropic. Right panels: inversion where anisotropy is imposed as an
unknown parameter at all depths. Modified from \citet{Alderetal2021}.}
\end{figure*}


Although the shear wave velocity is equally well resolved in both
cases, the anisotropic layer in the crust is better recovered with the
flexible scheme on the left. In the mantle, where the true model is
isotropic, the ensemble solution produced on the left panel includes a
large number of isotropic models, resulting in a narrower distribution
for radial anisotropy and a median value of the distribution equal to
unity. Conversely, in the case where anisotropy is imposed at all
depths (right panels), a wider distribution of anisotropy is obtained,
resulting in wider uncertainties. 

This comparison shows that in the case of a true isotropic model, a
flexible parameterization allows fitting data with simpler isotropic
models that are described by fewer parameters. The result of the
inversion is a distribution of models which is closer to the true
model.


\subsubsection{Exploiting the full wealth of information in posterior
distributions: example with azimuthal anisotropy}\label{sec3.2.2}

The previous example shows the benefits of using complex adaptive
parameterizations in Bayesian inversions. The solution is not a single
model, but an ensemble of models with variable\break parameterizations that
approximates a probability distribution defined in a multiple
dimensional space. As shown in Figure~\ref{fig8}, one simple way to
exploit this complex solution is to extract the mean or median value
for a given parameter at each depth. Standard deviations can also be
used, but more information is available about posterior covariances and
trade-offs, and inference about these quantities can be hard to make.
This wealth of information can be exploited through different angles of
analysis and \mbox{visualisation,} each one giving additional insight to
specific geophysical questions. 

We illustrate this in Figure~\ref{fig9} where Rayleigh wave dispersion
curves and their azimuthal variations measured from beamforming (see
Section~\ref{sec3.1.3}) are inverted at depth following the method of
\citet{Bodinetal2016}. Here, the number of layers is variable and each
layer is either isotropic and described solely by its isotropic
shear-wave velocity, or azimuthally anisotropic and described by three
parameters: (1)~the isotropic shear-wave velocity; (2)~the peak to peak
level of azimuthal anisotropy, and (3)~the direction of the horizontal
fast axis relative to the north \citep{RomanowiczYuan2012}.


\begin{figure*}
\vspace*{-3pt}
\includegraphics{fig09}
\vspace*{-2pt}
\caption{\label{fig9}Depth inversion of Rayleigh wave surface wave
dispersion curves with azimuthal variations for a subarray located in
the Dinarides. The solution is a large ensemble of models with
different parameterizations (different numbers of isotropic and
anisotropic layers). (a)~Probability distribution of the amplitude. 
(b)~Probability distribution of the direction of the fast axis of
anisotropy. The black line in (a) shows the average amplitude, and the
horizontal black lines in (a) and (b) indicate the limits of the three
layers over which we integrate anisotropy. (c)~Density plot of
anisotropy parameters integrated over the three different depth ranges.
The radius in the polar plots indicate amplitude and the angle gives
direction of the fast axis of anisotropy. The colour scale indicates
the normalised probability of a given combination of amplitude and fast
direction. Panels (d) and (e) show the marginal distributions of
${\delta}V_{\mathrm{s}}$ (d) and of $\theta_{2}$ (e) for each
of the three layers. The vertical scale is chosen such that the surface
of the blue area is 1, that is the fraction of models with a given
amplitude or ($\theta_{2}$) range corresponds to the
blue area within that range. Modified from \citet{Soergeletal2023}.}
\vspace*{-3pt}
\end{figure*}


Figure~\ref{fig9} shows an example of inversion of a set of beamforming
measurements below the Dinarides (location:
${\sim}$45.3\textdegree{}N--16.3\textdegree{}E). The ensemble solution is
shown in the left panels, where the probability of a given amplitude
range and direction range for shear-wave anisotropy is shown at each
depth. However, this way of displaying the posterior distribution does
not show correlations and trade-offs between anisotropy amplitude and
direction. Such plots do not show whether the sampled profiles contain
thin isotropic layers, thin layers with strong anisotropy, or thick
layers with small anisotropy.

To extract some key features from this complex probabilistic solution,
\citet{Soergeletal2023} proposed to show the distribution of integrated
anisotropy over a given depth range. The integration was defined as in
\citet{RomanowiczYuan2012}, and carried out over three specific depth
ranges with a meaningful sense (upper crust, lower crust, uppermost
mantle). The distribution for the amplitude and direction of the
anisotropy integrated over the three layers is shown in the right
panels. 

This visualisation makes it possible to see the correlations between
fast direction and amplitude of anisotropy in each depth interval. Note
that isotropic layers decrease the average level of anisotropy.
Additionally, integrating anisotropy over several layers with different
fast directions can also lead to anisotropy amplitudes smaller than
expected. This explains the differences between the distribution of
local anisotropy (left panels), and integrated anisotropy (right
panels). 


\section{New information on the greater Alpine\newline region}\label{sec4}

One of the main motivations behind the AlpArray initiative was to bring
the resolution of seismic tomography closer to the spatial density of
geological data over the entire Alpine belt, with the aim of increasing
resolution of three-dimensional geological models at lithospheric
scale. The methodological developments described in Section~\ref{sec2}
have made a \mbox{valuable} contribution to this aim, in particular by
providing new insights into the geometry of the crust-mantle boundary.
A team of geologists and geophysicists is now working on exploiting the
shear-wave velocity models derived by \citet{Nouibatetal2022a} and
\citet{Nouibatetal2023} to build a 3-D structural model of the Western
Alps \citep{Baderetal2023}. All available geometrical (digital
elevation model), geological (structural and geological maps) and
geophysical (3-D $V_{\mathrm{p}}$ and $V_{\mathrm{s}}$ models, DSS profiles, Moho models from
seismic and gravity modeling/inversion, etc.) have been integrated and
mixed in a geomodeling software, which provides a common framework and
checks the geometrical coherence of geological\break interpretations of
geophysical data \citep{Calcagnoetal2008}. The resolution of
\citet{Nouibatetal2022a}'s model was sufficiently fine to allow
\citet{Sonnetetal2023} to interpret a lateral $V_{\mathrm{s}}$ change in the
subducted European lower crust as indication for the transition from
amphibolite to granulite based on petrophysical data on Alpine rocks.

A depth map of the iso-velocity surface $V_{\mathrm{s}}= 4.2$~km/s extracted from
the model by \citet{Nouibatetal2023} is shown in Figure~\ref{fig10}. As
explained in Section~\ref{sec2.4}, a composite depth map of velocity
surfaces 3.8, 4.1 and 4.3 km/s would be a better proxy of the three
Moho boundaries for the Ivrea body, the \mbox{Ligurian} basin, and Eurasia and
Adria outside the Ivrea body region. The $V_{\mathrm{s}}= 4.2$ km/s surface is
therefore only an easy-to-calculate compromise. Since the velocity
gradient is generally strong between 4.1 and 4.3 km/s,
Figure~\ref{fig10} well illustrates lateral changes in Moho depth
beneath the Alps revealed by our ambient noise tomography studies. The
purple arrows highlight an ${\sim}$8 km step in the European Moho with
SSW--NNE orientation beneath the external crystalline massifs of the
Western Alps, from Pelvoux to Mont Blanc. This orientation suggests
that the Moho step might be a major lithospheric structure inherited
from the Variscan orogeny. It had never been imaged before, even by the
ECORS-CROP deep reflection profile, which showed reflections from the
lower crust only west of the Belledonne external crystalline massif
{[}e.g., \citealp{Nicolasetal1990}; see also Supplementary Figure~S3 in
\citet{Pauletal2022}{]}. The $V_{\mathrm{s}}$ models of  
\citet{Zhaoetal2020,Nouibatetal2022a}, 
and  \citet{Nouibatetal2023} have confirmed
the subduction of the European lithosphere beneath Adria, which is
documented in Figure~\ref{fig10} by the dark blue areas (depth
$>$ 60~km). The CIFALPS cross-sections in Figures~\ref{fig3}d--f
document this continental subduction more precisely, with a European
Moho reaching 75--80~km depth. Figure~\ref{fig10} also highlights the
strong and rapid Moho depth changes along the belt strike, with maximum
depths located right beneath the Insubric line (IL in
Figure~\ref{fig10}, also named Periadriatic fault further east) which
marks the western and northern\break boundary of the undeformed Adria
microplate \citep{Handyetal2010}. The 3-D detailed geometry and
internal structure of the Ivrea body are also important new findings of
our ambient noise tomography studies. For example, the ``inverted Moho''
highlighted by receiver function sections (InvM in
Figure~\ref{fig3}b, dotted thick grey line in Figure~\ref{fig3}e) is
explained by the overthrusting of the Ivrea mantle slice of Adriatic
origin onto the lower velocity European crust. This boundary was imaged
a few tens of km north of the CIFALPS line by wide-angle reflections of
the ECORS-CROP complementary experiment
\citep{ECORS-CROPDeepSeismicSoundingGroup1989}. It was interpreted as
the top of a second mantle slice beneath the Ivrea body, while the
polarity of converted signals in receiver function studies and ambient
noise tomography studies have proven that this reflector is the base of
the Ivrea body (single) mantle slice \citep{Zhaoetal2015,Pauletal2022}.


\begin{figure*}
\vspace*{5pt}
\includegraphics[scale=0.95]{fig10}
\vspace*{2pt}
\caption{\label{fig10}Depth map of the iso-velocity 
surface $V_{\mathrm{s}}=4.2$
km/s, extracted from the wave-equation ambient noise tomography model
of \citet{Nouibatetal2023}. The continuous and dotted black lines are
the main geological boundaries and faults, similar to
Figure~\ref{fig1}. The purple arrows highlight the Moho step located
beneath the external crystalline massifs of the Western Alps. The
dashed white line outlines the top of the high-velocity, peridotitic
core of the Ivrea body, which appears as a dark blue spot in
\mbox{Figures~\ref{fig3}d--f.} Ad: Adria; CA: Central Alps; EA: Eastern
Alps; Eu: Eurasia; IL: Insubric line; Li: Ligurian basin; WA: Western
Alps. Modified from \citet{Nouibatetal2023}.}
\end{figure*}


In the Ligurian basin, the wave-equation ANT model with
elastic--acoustic coupling by \citet{Nouibatetal2023} has provided new
images of the basin margins and new estimates of sediment thickness
with potential inferences on tectonic models of basin formation. On
Ligurian Moho depth, the recently published ANT $V_{\mathrm{s}}$ model by
\citet{Magrinietal2022} exhibits significant differences 
from the $V_{\mathrm{s}}$
model by \citet{Nouibatetal2022b}, with 15--18~km crustal thickness
along the basin axis in \citet{Magrinietal2022}'s model and 12~km in
\citet{Nouibatetal2023}'s. As explained in Section~\ref{sec2.4}, the
thinner crust of \citet{Nouibatetal2022b} is fully consistent with the
P-wave velocity model derived from the refraction-reflection profile
LOBSTER-P02 of \citet{Dannowskietal2020} along the Ligurian basin axis.
The difference in Moho depth estimates between
\citet{Nouibatetal2022b,Nouibatetal2023} and \citet{Magrinietal2022}
may be due to different ray coverages in the Ligurian basin, as
\citet{Magrinietal2022} did not use AlpArray OBS data. 

Very little information was available on the anisotropy of the crust
and upper mantle beneath the Alpine region prior to the works described
in Section~\ref{sec3}. \citet{Kastleetal2022} and
\citet{Soergeletal2023} have mapped azimuthal anisotropy of Rayleigh
wave phase velocity using ambient noise records of the AASN and
permanent seismic stations, and two different methods, Eikonal and
beamforming. Their azimuthal anisotropy maps are broadly similar at all
periods, except at the periphery of the study area, probably due to
differences in ray coverage\break (see examples at 15 and 40~s period in
Figures~\ref{fig5} and~\ref{fig7}). The Bayesian inversions for
azimuthal anisotropy distribution with depth by \citet{Soergeletal2023}
have shown that the anisotropic structure cannot be easily inferred
from maps at individual periods. Spatially coherent anisotropy patterns
are only visible in the upper half of the crust, with fast-velocity
directions mostly parallel to the strike of the belt and amplitudes of
1--2\% \citep[see Figure~14 in][]{Soergeletal2023}. By contrast,
fast-velocity directions are mostly perpendicular to the belt in the
lower half of the crust and uppermost mantle, with strongly varying
amplitudes and no large-scale spatial pattern. This contrast between
the shallow and deep layers suggests that Alpine deformation has 
only\break impacted the upper crust, through oriented 
crack and fractures, while
the lower crust and upper mantle bear the imprint of more ancient
processes. A single region located northwest of the Jura mountains
displays coherent NE--SW fast velocity directions from the upper crust
to the upper mantle, and rather strong amplitudes (${\sim}$1\%) in the
lower crust. This orientation suggests that the observed anisotropy in
that area outside the Alps may be of Variscan origin, and unaffected by
the Alpine orogeny. \citet{Soergeletal2023} highlight the general
disagreement between the fast-velocity directions they measured in the
upper mantle and the fast directions measured from the splitting of
core-refacted XKS teleseismic phases \citep[e.g.,][]{Heinetal2021}. This
may indicate that the source of XKS splitting is located deeper than
the lithospheric mantle, either in the asthenosphere or in the
subduction slabs.
\pagebreak

When applied to Rayleigh and Love dispersion curves measured in the
Alpine region, the flexible parameterization of Bayesian inversion
allowed \citet{Alderetal2021} to produce maps of radial anisotropy
differing from previous large-scale studies that suggested the presence
of significant radial anisotropy everywhere in the European crust and
shallow upper mantle \citep[e.g.,][]{Zhuetal2015}. Instead, they observed
that radial anisotropy is mostly localised beneath the Apennines while
most of the crust and shallow upper mantle is isotropic in other parts
of Europe. Thanks to synthetic tests, they attributed this difference
to trade-offs between radial anisotropy and thin (hectometric) layering
in previous studies based on least-squares inversions and long period
data (${>}$30~s). In contrast, \citet{Alderetal2021}'s approach involved
a massive dataset of short period measurements and a Bayesian inversion
that accounts for thin layering. They showed that the positive radial
anisotropy ($V_{\mathrm{SH}} > V_{\mathrm{SV}}$) observed in the lower crust of the Apennines
could not result from thin layering, but rather from ductile horizontal
flow in response to the strong flexure of the Adriatic plate induced by
doubly-vergent subduction. 

\section{Conclusions}\label{sec5}

The Alpine broadband seismic networks, including permanent stations and
the temporary AlpArray seismic network have provided an optimal dataset
for the development, improvement, and application of
ambient-noise-based imaging methods. High station density and
homogeneous coverage were key elements for these methodological
improvements. 

The overarching goals of our work were to improve observations when
needed, move towards a probabilistic framework for inversion and
interpretation, and adapting full waveform modelling in the forward
problem. On the observational side, the work focused on (a)~developing
suitable techniques to combine OBS and land observations through higher
order correlations, (b)~improving our understanding of systematic errors
made when measuring surface wave azimuthal anisotropy, and (c)~applying
and adapting beamforming techniques to seismic noise correlations to
estimate surface wave anisotropy. 

For the inversions, a Bayesian framework was applied both to 2-D
inversions for group velocity maps \citep{Nouibatetal2022a}, and to 1-D
depth inversions for elastic structure
\citep{Luetal2018,Zhaoetal2020,Alderetal2021,Soergeletal2023}.
Particular efforts \citep{Alderetal2021,Soergeletal2023} were made on
the extraction of robust anisotropic\break parameters which are meaningful in
a geologic/{\ubreak}tectonic sense. Wave-equation ambient noise tomography
studies by \citet{Luetal2020} and \citet{Nouibatetal2023} have further
improved ANT models by accounting for the physics of seismic wave
\mbox{propagation.} The introduction of elastic--acoustic coupling at the
seabed has enabled \citet{Nouibatetal2023} to improve the imaging of
shallow layers below the sea. All of these developments led to improved
models of the crust and uppermost mantle on the scale of the greater
Alpine area (see Section~\ref{sec4}).

The link between the overarching goals is naturally increasing because
of the strong link between observational errors and the estimation of
model resolution within a probabilistic inversion framework. In such a
framework, the inversion provides not a single Earth model, but an
ensemble of models obtained whilst taking into account uncertainties of
inverted data. In this context, reliable error estimates of observables
become objectively tied to the estimation of model uncertainties
\citep{Gallagheretal2009,Sambridgeetal2013}. Also, the ensemble of
models can be explored based on scientific questions (for example ``What
is the probability that the lower crust beneath the Apennines has
significant radial anisotropy of crystallographic origin?'')
\citep{Curtis1999,Meieretal2007,ZhangCurtis2021}. 

Note also that the Bayesian framework allows to include independent
information from geology, mineralogy, or other geophysical methods such
as gravimetry. This can be done through the definition of the a priori
distribution, or through the joint inversion of different datasets.
Including other data types sensitive to P-wave structure (body waves)
would be very useful for geological interpretation, as P and S waves
have different sensitivities to different rock types, for example to
hydrated minerals. In this context, full waveform inversion has the
benefit to model the entire wavefield without having to preliminary
extract and separate different data types (arrival times, phase
velocities, splitting parameters, etc.), and three-dimensional
variations of the entire\break elastic tensor can be directly reconstructed,
including $V_{\mathrm{p}}$, $V_{\mathrm{s}}$ and anisotropy. However, full waveform tomography in a
fully non-linear Bayesian framework is still a prospect for the future,
due to the computational cost of the forward model.

Within the framework of the AlpArray experiment, other noise-based
methods were developed. For crustal imaging, it is becoming feasible to
image coda-Q \citep{Soergeletal2020}, overcoming the difficulty of
obtaining reliable amplitudes of noise correlations. This method makes
it possible to estimate coda-Q at longer periods than those accessible
through local earthquake estimates, but the method needs to be further
developed to fully understand the meaning and reliability of the
observed coda-Q values. The dense data coverage of the Alps also served
the purpose of improving methods to extract body waves (reflections
from mantle discontinuities) from seismic noise correlations
\citep{Pedersenetal2023,Luetal2023}. 


The methods developed within the framework of the AlpArray experiment
can naturally be applied and further improved within the framework of
other dense and large-scale seismic arrays, such as the recently
started AdriaArray experiment which significantly extends the area of
high station density around the Adriatic plate, from the French Massif
Central in the west to the Carpathians in the east
[\citealp{Kolinskyetal2023};
\url{https://orfeus.readthedocs.io/en/latest/adria_array_main.html}{]}.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds
from any organization that could benefit from this article, and have
declared no affiliations other than their research organizations.

\section*{Funding}
A large part of the studies described here were conducted
within the AlpArray-FR project funded by Agence Nationale de la
Recherche, France (contract ANR-15-CE31-0015), and Labex OSUG@2020
(Investissement d'Avenir, ANR-10-LABX-56). TB, CA and DS were funded
by the European Research Council under the European Union Horizon 2020
research and innovation program (grant agreement 716542 --
TRANSCALE). AP's work on the Alps is currently funded by project
LisAlps (contract ANR-20-CE49-0007). EK is funded by the German Science
Foundation DFG (SPP-2017, Project Ha 2403/21-1).

\section*{Acknowledgments}
We are grateful to our colleagues, Laurent Stehly, Romain Brossier,
and Eric Debayle who co-supervised YL's, CA's and AN's Ph-D theses with
AP and TB, and have made major contributions to much of the work
presented here. Liang Zhao (IGG-CAS Beijing) is also acknowledged for
his decisive contributions to this work. This paper was partly
conceived during the Carg\`{e}se 2022 school on ``Passive imaging and
monitoring in wave physics: from seismology to ultrasound'', in which
AP took part. We would therefore like to thank the organizers, in
particular Michel Campillo and Nikola\"{i} Shapiro for encouraging us
to write it. We acknowledge detailed and constructive comments by two
anonymous reviewers that helped improving the manuscript, and editorial
handling by Andrew Curtis. Without high-quality, open data, such work
is impossible. We therefore gratefully thank the operators of the
European permanent seismic networks who make their data available
through EIDA (\url{http://www.orfeus-eu.org/eida}), the AlpArray
Seismic Network Team who operated the AASN (Z3-2015 network) and the
CIFALPS Team who operated the CIFALPS experiments (YP-2012 and XT-2018
networks). 


%\back{}

%\refinput{crgeos20230948-reference.tex}
\bibliographystyle{crgeos}
\bibliography{crgeos20230948}

\end{document}