Nishida et al.
Progress in Earth and Planetary Science(2024) 11:4
https://doi.org/10.1186/s40645-023-00605-8

Progress in Earth andPlanetary Science

Open Access

REVIEW

Ambient noise multimode surface wave
tomography
Kiwamu Nishida1* , Ryota Takagi2 and Akiko Takeo1

Abstract
Seismic techniques using earthquakes are powerful tools for exploring the Earth’s internal structure. However,
the earthquake distribution limits the spatial resolution. In recent years, ambient noise surface wave tomographyusing ambient seismic wave field has resolved these limitations. A typical ambient seismic wave field is microseisms
excited by ocean swell activities. Ambient noise surface wave tomography is a technique in seismic interferometrythat extracts seismic wave propagation between pairs of stations by cross-correlating the seismic records. The crosscorrelation function can be interpreted as an impulsive response at a station with a virtual source at the other station.
This technique became standard with the development of modern dense seismic networks. This paper reviews a theory of seismic interferometry for ambient noise surface wave tomography and procedures for practical data processing to calculate cross-correlation functions. The tomographic method typically consists of four steps: (1) the construction of reference 1-D models, (2) phase velocity measurements for each path, (3) 2-D phase velocity inversions, and (4)the construction of a 3-D S-wave tomographic model obtained from series of local 1-D inversions at all the grids. This
paper presents the feasibility of multimode surface wave dispersion measurements for improving depth resolution.
Keywords Ambient seismic noise, Multimode surface wave, Seismic interferometry, Surface wave tomography

1 Introduction
Seismic techniques are powerful tools for exploring the
Earth’s internal structure from local to global scales.
Earthquakes are a primary illumination source of the
Earth’s interior. Seismic tomography using earthquakes
revealed the lateral heterogeneities of the Earth on various scales (e.g., Romanowicz 2003; Thurber and Ritsema2015). Although the method has produced many significant results, it has inherent limitations attributed to theearthquake distribution. Earthquakes occur only in tectonically active areas, and large earthquakes do not occur*Correspondence:
Kiwamu Nishidaknishida@eri.u-tokyo.ac.jp
1
Earthquake Research Institute, The University of Tokyo, 1‑1‑1 Yayoi,
Bunkyo‑ku, Tokyo 113‑0032, Japan
2
Research Center for Prediction of Earthquakes and Volcanic Eruptions,
Graduate School of Science, Tohoku University, 6‑6 Aza‑Aoba, Aramaki,
Aoba‑ku, Sendai 980‑8578, Japanfrequently. The source and station distributions limit the
spatial resolution of the tomographic images.
Even on seismically quiet days, the Earth persistently
oscillates because ocean swells excite seismic waves. They
are known as microseisms (e.g., Nishida 2017), which is akind of ambient seismic wave field (also known as ambient seismic noise). Seismic interferometry (SI) provided
a clue to overcome the limitations due to the earthquake
distribution because SI produces a virtual earthquake
record from the ambient seismic wave field. An application of SI for imaging Earth’s internal structure is known
as ambient noise surface wave tomography (ANT). Thissection briefly introduces SI, ambient seismic wave field,
and ANT followed by the scope of this paper.
1.1 Brief history of seismic interferometry

In seismology, an origin of SI is Aki’s Ph.D. thesis (Aki1957), inspired by the seminal book Cybernetics (Wiener 1947). He proposed the SPatial AutoCorrelation(SPAC) method; later studies show that the method is© The Author(s) 2024. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://creativecommons.org/licenses/by/4.0/.

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

mathematically equivalent to SI (e.g., Chávez-Garía andLuzón 2005; Harmon et al. 2010; Tsai and Moschetti2010) under certain conditions. Aki’s idea had not been
paid attention until the 1980s. After 26 years, a group at
Hokkaido University (Okada and Sakajiri 1983) developed a survey method for shallow structure using microtremor, now known as the microtremor survey. After thiswork, this method became popular for surveying subsurface structures in seismic engineering (e.g., Cho et al.
2006), which was also extended to multimode inversion(Tokimatsu et al. 1992).
A key for SI is a cross-correlation function (CCF)
between a pair of stations. The CCF exhibits an impulsive response at a station by a virtual source at the otherstation. The idea of SI dates back to the 1950s in different research fields: acoustic (Eckart 1953), ocean acoustic (Cox 1973), seismic exploration (Claerbout 1968),
and seismology (Aki 1957). Although the ideas were proposed independently, they are mathematically identicalunder certain conditions (e.g., Chávez-Garía and Luzón
2005; Harmon et al. 2010; Tsai and Moschetti 2010). SI
is applied not only to the Earth but also to experimental studies, such as ultrasonic waves (Lobkis and Weaver
2001), building responses (Snieder and Wapenaar 2010;
Nakata et al. 2013), helioseismology (Gizon and Birch
2002; Duvall et al. 1993; Hanasoge et al. 2016) and ocean
acoustic (Cox 1973). Later studies (e.g., Godin 2007;
Snieder and Larose 2013) show that the principle of SI
is akin to the fluctuation–dissipation theorem in physics (e.g., Callen and Welton 1951), which was generalized
as linear response theory by Kubo (1957). The theorem
establishes a connection between the response to an
external force and the CCF of specific fluctuating properties when the system reaches thermal equilibrium.
1.2 Ambient seismic wave field

Ocean wave activities excite ambient seismic wave fields,
even on seismically quiet days. Based on the types ofocean surface gravity waves and excitation mechanisms,
we classified ambient seismic wave fields into (1) seismic
hum from 1 × 10−3 to 0.02 Hz, (2) primary microseisms
from 0.02 to 0.1 Hz, and (3) secondary microseisms
between 0.1 and 0.5 Hz (Nishida 2017).
The ocean surface gravity waves can be classified as
ocean infragravity (IG) waves below 0.02 Hz and oceanswell above the frequency. IG wave is a shallow-water
wave, whereas the ocean swell is physically a deep-water
wave. The pressure fluctuations of the IG waves reachthe seafloor in the pelagic and coastal regions, whereas
those of the ocean swells cannot reach the seafloor in the
pelagic regions.
IG waves excite background Love waves and Rayleigh waves predominantly from 1 to 20 mHz, known asPage 2 of 42

seismic hum, through the topographic coupling on theseafloor (Nishida et al. 2008b; Fukao et al. 2010; Nishida2013). Ocean swell activities excite microseisms from0.05 to 0.5 Hz (Nishida 2017). Primary microseisms are
excited by shoaling waves through the topographic coupling on the sea floor (e.g., Hasselmann 1963; Ardhuin
et al. 2015). Secondary microseisms are excited by oceanswell activities in both the pelagic and coastal regions
through the nonlinear interaction (Longuet-Higgins
1950; Hasselmann 1963). In the history of earthquake
seismology, secondary microseisms have been a prominent noise source for earthquake signals. Historicalseismological studies took different approaches depending on two frequency bands below and above the typicalfrequency of secondary microseisms. This frequency isalso important in characterizing the physical nature ofseismic waves (Aki 2003). Above the frequency, strong
lateral heterogeneities in the crust and sediment makethe seismic wave field complex on a regional or global
scale (see Sect. 6.2 for details), stimulating the development of a stochastic approach. Below the frequency, thewaveform can be reproduced in a deterministic manner.
SI turns stochastic wave fields into a deterministic signal
as a virtual seismic record, which can be applied to seismic tomography. Thus, the operation of cross-correlating seismic wave fields, or SI, bridges the stochastic anddeterministic approaches in seismology, which had been
incompatible, as noted by Aki (2003).
1.3 ANTCampillo and Paul (2003) demonstrated that SI can
extract virtual seismic records from the ambient seismic
wave field. This paper describes the analysis of the coda
wave, a diffusive wave field, in Mexico. Cross-correlating
the seismic records of coda waves extracted clear Love
and Rayleigh wave propagations between every pair ofstations. Although this paper analyzed earthquake coda,
the method can be applied to ambient seismic wave
fields. Shapiro and Campillo (2004) subsequently showedthat the CCFs of the microseisms also exhibited clear
surface propagations between every pair of stations. This
result suggested the possibility of seismic imaging without earthquakes, now known as ANT.
Shapiro et al. (2005) achieved the milestone of the
ANT studies. They demonstrated the feasibility of ANT
using a modern dense seismic network. After this study,
ANT became a standard technique for dense broadband observations because ANT does not require waiting for earthquakes, i.e., a long observation period. Afterthe paper, ANT was applied to many regions: the USA(e.g., Bensen et al. 2007; Moschetti et al. 2007; Liang andLangston 2008; Lin et al. 2008; Yang et al. 2008), Australia(e.g., Saygin and Kennett 2010), Europe (e.g., Yang et al.

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

2007), China (e.g., Zheng et al. 2008), Japan (e.g., Nimiya
et al. 2020; Nishida et al. 2009), and the global scale(Nishida et al. 2009; Haned et al. 2016). With the increase
in seismic stations with continuous observations at many
seismic stations in the 2000s, ANT has become a standard analysis method. Today, ANT is first applied after
deploying a dense seismic array.
1.4 What is covered/not covered in this review

There are already many review papers on SI and
ANT (e.g., Snieder and Larose 2013; Wapenaar et al.
2010a, 2010b) and textbooks (e.g., Schuster 2009, Ritzwoller and Feng 2019, Nakata et al. 2019, and Chapter 10
in Sato et al. 2012). This review focuses on a consistenttheoretical treatment and systematic comparison among
different phase velocity measurement methods.
This review also focuses on multimode measurements. Recent developments in dense arrays enable us
to extract multimode dispersion (e.g., Spica et al. 2018;
Chmiel et al. 2019; Savage et al. 2013; Jiang and Denolle
2022; Socco et al. 2010). Multimode inversions have significantly improved depth resolution, although the dispersion of fundamental mode branches alone has poor
vertical resolution. As multimode dispersion measurement has inherent difficulties, this review aims to provide
clues for practical applications.
Because SI does not require earthquakes, it accelerates the progress in monitoring the temporal change inseismic velocities. Monitoring the temporal change with
events requires repeating earthquakes or repeating active
sources. In most cases, they are not realistic. In the lastten years, SI has become a standard technique for monitoring the temporal change in seismic velocities associated with environmental origins (e.g., precipitation,
groundwater, and temperature), volcanic eruptions, andearthquakes (e.g., Sens-Schönfelder and Wegler 2006a;
Wegler et al. 2007; Brenguier et al. 2008a, b, Wang et al.
2017). Because this review does not cover such topics, seereview articles (e.g., Sens-Schönfelder and Wegler 2011;
Obermann and Hillers 2019) for further information.
This paper does not cover conventional surface wave
tomography using earthquake data. Because many excellent reviews are already available (e.g., Romanowicz
2003, 2020, 2021; Laske and Widmer-Schnidrig 2015;
Levshin et al. 2018; Barmin et al. 2001), please refer tothose references.
This paper also does not cover attenuation tomography using CCFs of ambient seismic noise. Some studies
inferred the attenuation structure from the amplitude
information of CCFs (Prieto et al. 2009, 2011; Lin et al.
2011). Although the effects of source heterogeneities(e.g., Tsai 2011) can bias the estimation, recent developments to extract amplitude information from ambient

Page 3 of 42

noise cross-correlations (e.g., Lin et al. 2012b; Zhou et al.
2020; Liu et al. 2021) have made it reliable. Because this
topic is beyond the scope of this paper, we refer only tothe aforementioned papers.
Section 2 summarizes a theory of SI for ANT. Section 3 describes the procedures for practical data processing to calculate CCFs. Ambient noise multimodesurface wave tomography (multimode ANT) consists offour steps (Fig. 1). The first step (Sect. 4) is the measurement of multimode surface wave dispersion for a local
seismic array. The dispersion curves are also crucial for
constructing local 1-D structures, which can be an initialmodel of ANT. The second step (Sect. 5) is the dispersion
measurements for each path. The third and fourth steps(Sects. 6 and 7) are how to infer 3-D seismic velocity
models from the dispersion measurements. The procedures consist of two steps: (1) the 2-D inversion of phase/group velocities in Sect. 6 and (2) a local 1-D inversion
at each grid in Sect. 7. Figure 1 shows such procedures
based on our previous studies (Nishida et al. 2008a;
Nagaoka et al. 2012; Takagi and Nishida 2022; Takeo
et al. 2022; Yamaya et al. 2021). This review will focus onand compare our previous studies with other studies.

2 A brief review of SI
Theories of SI originate from various backgrounds andhave developed independently. Although they requiredifferent assumptions in different settings (e.g., an open
or closed system), they are closely related. Theoretically,
it is natural to consider a closed system for the global
scale, whereas it is natural to consider an open system forthe regional or local scale. This section provides an overview of the theories of SI in a closed system and an open
system to understand the physical pictures. To model
CCFs for ANT, we formulate synthetic CCFs for the
homogeneous source distribution in a 2-D homogeneous
medium. However, the heterogeneous source distribution realistically causes an apparent travel-time anomaly
from the synthetic CCF for the homogeneous source distribution. We evaluate such biases based on the analytic
formulation.
2.1 SI in a closed system

Here, we consider SI in a closed system. In the case of afinite body, we evaluate CCFs based on a normal mode
approach (Lobkis and Weaver 2001). Because the Earthis a finite-size sphere, this approach is also feasible formulti-orbit propagations on a global scale (Nishida et al.
2002, 2009). For simplicity, this subsection describes ascalar 1-D case, but it can be easily extended to elastic
2-D and 3-D cases.

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Page 4 of 42(a) 1st step: secion 4
Construction of a reference
dispersion curves of an array(b) 2nd step: section 5
Phase velosity anomaly
for station pairs

LagtimeLagtime

Reference dispersion curves

Dispersion curves for each path
Phase vel.

Phase vel.
Frequency(c) 3nd step: section 6
Phase velosity tomography
0th RayleighFrequency

Frequency

Love wave

Phase vel.

Rayleigh wave

Love wave

Phase vel.

Rayleigh wave

ZZ, ZR,
RR, TT

Distance

Distance

ZZ, ZR,
RR, TT

Frequency(d) 4th step: section 7
1D inversion at all the grid points

1st Rayleigh3D Vs structure
1D inversion

Fig. 1 a Construction of reference 1-D models. b Phase/group velocity measurements for each path. c 2-D phase/group velocity inversions
for multimodes. d A 3-D S-wave tomographic model is obtained from a collection of local 1-D inversions at all the grids. Modified from Takagi
and Nishida (2022)

2.1.1 Repeating seismic experiments

Virtually, we consider repeating seismic experiments ina closed system of a perfect elastic body. Before the initialtime t = 0, the body did not deform, and a random force
f k (x) was applied to the body at t = 0 in the kth experiment, where x is the spatial location. We consider a finite
body from x = 0 to x = L and impose rigid or free boundary conditions at both ends. The equation of motions isgiven byρ(x)∂ 2 uk (x, t)∂ 2 uk (x, t)=
κ(x)+ f k (x)δ(t),
∂t 2
∂x2

(1)

where κ is the elastic constant, ρ is the density and uk isthe displacement at the kth experiment. After the force
is applied, the displacement is measured at the stations.
These experiments are repeated K times. Note that wecannot consider a persistent force in this system because
no attenuation leads to an infinite increase in the amplitude over time.
We evaluated the displacement uk (x, t) by a convolution between the force f k (x) and the Green’s function.
The Green’s function can be written in terms of normal
mode theory (Dahlen and Tromp 1998) as

Nishida et al. Progress in Earth and Planetary Science

G(x, x′ ; t) =

 Un (x)Un (x′ )nωn(2024) 11:4

sin(ωn t), t ≥ 0,

Page 5 of 42(2)

where ωn is nth eigenfrequency and Un is nth eigenfunction, which satisfies orthonormality:
 L
ρ(x)Un (x)Un′ (x)dx = δnn′ .
(3)0Then, the displacement can be represented by the
convolution

uk (x, t) =



Akn Un (x)nsin(ωn t)
,
ωn

(4)

L

where Akn = 0 Un (x)f k (x)dx . We consider the energy
partition of the modes for a random external force in thenext section.
2.1.2 Energy partition of modal energy

This subsection discusses the energy balance of each
mode: how the work done by external forces is distributed to the kinematic and elastic energy. First, we evaluate the work done by the external force can be given by
 ∞ Lvnk (x, t)f k (x)δ(t)dxdt = (Ank )2 ,
(5)
0

0

where the particle velocity vnk of nth mode in the kth
experiment is written by vnk = Akn Un (x) cos(ωn t).
The kinetic energy of nth mode Tn can be evaluated byintegrating the kinetic energy density in space as
 L
1
ρ k 2
(vn ) dx = (Ank )2 cos2 (ωn t).
Tn =
(6)
2
0 2
The elastic energy Vn, on the other hand, can be evaluated by integrating the strain energy density in space. The
partial integral with the boundary condition leads to
2
 L 
1
κ
k ∂ Un sin(ωn t)
An
dx = (Ank )2 sin2 (ωn t).
Vn =
∂x
ωn
2
0 2
(7)
The total energy Tn + Vn isTn + Vn = (Akn )2 ,(8)

which balances the work done by external forces.
Here, we consider an excitation by random force Fik at
location xi as

f k (x) =

I−1

i=0

Fik δ(xi ),(9)where I is the number of the force. We also assume that
Fik (x) is white noise as�Fik Fik′ �k = F̄ 2 δii′ ,(10)where k is ensemble average with respect to k. To simplify the problem, we consider the constant density ρ0.
The expected value of the cross-correlation of An can beevaluated as follows.
I F̄ 2
I F̄ 2 Lk k
Un (x)Un′ (x)dx =δnn′ ≡ E δnn′ ,An An′ ∼
k
Lρ0 0
Lρ0
(11)
where E is the modal energy. The amplitudes of the different modes An do not correlate with each other, andthe total energy of each mode is distributed equally. The
expected value of the modal energy is constant for eachmode when the external force f is white noise. When asystem meets this condition, we call the state the equipartition of energy. We note that the fluctuation–dissipation theorem in physics (Callen and Welton 1951)
requires a thermal equilibrium that satisfies the equipartition of energy. Equation 10 can be derived from theprinciple of equal a priori probabilities in the case.
2.1.3 CCFs under the equipartition of energy

Here we define a CCF φ k (x1 , x2 ; τ ) between uk (x1 ) anduk (x2 ) under the equipartition of energy as,

1 T k
k
φ (x1 , x2 ; τ ) ≡ lim
u (x1 , t)uk (x2 , t + τ )dt,
T →∞ T 0(12)
where x1 shows the location of station 1, and x2 showsthat of station 2. We note that there are two types of CCF
definitions. The sign of the second term of the other typeis flipped.
The ensemble average of φ(x1 , x2 ; t) over
 K time experiments is defined by the ensemble average φ k (x1 , x2 ; τ ) k .
Then, the CCF can be represented bydE
φ(x1 , x2 ; τ ) = − (G(x1 , x2 ; τ ) − G(x2 , x1 ; −τ )),
dτ
2
(13)
which relates the CCF to Green’s function (Snieder
2004). The frequency dependence of this equation differs
slightly from Aki’s formulation (e.g., Haney et al. 2012).
Although the formulation showed that the Hilbert transform of the CCF can be related to Green’s function, thedifference can be attributed to spectral normalization.
Compared to CCFs for real data, the biggest problem
is the assumption of the equipartition of energy. Because
the excitation sources are distributed near the surface,fundamental modes dominate the observed wave field.

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Energy is not partitioned equally among different mode
branches. Indeed, the observed dominance of fundamental modes (e.g., Nishida 2017) shows that the energy is
not equally distributed in the radial direction. When considering a single-mode branch, the energy is distributed
equally in the horizontal direction if the external force iswhite noise distributed on the entire surface.
Although the theory assumed K experiments, practically, only one observation is possible. Therefore, theensemble averages need to be replaced by time averages
to calculate the CCFs of observed data for persistentexternal forces (e.g., ocean waves). However, no attenuation in a closed system causes the problem of diverging
amplitudes without attenuation; that is, the seismic wave
field in a closed system never meets the equilibrium for
persistent sources. Physically, it is natural to consider theequilibrium between the energy dissipation due to attenuation and the work performed by external forces (Kobayashi and Nishida 1998; Fukao et al. 2002).
We note that the term Green’s function is often used
in a less mathematically rigorous sense in studies on SI.
Mathematically, a CCF converges to the Green’s function
only in limited cases. A similar situation occurs in quantum field theory. The Green’s function refers to CCFs,
although it does not satisfy the mathematical definition(Zagoskin 2014).

Page 6 of 42

open system. When external force acts on the body at alocation r and time t = 0, the 2-D wave equation is givenby∞∂ 2 ψ(r, t)= C(r)∇ 2 ψ(r, t) +δ(ri − r)δ(t)fi
2
∂t
i=0

(14)
where C(r) is phase velocity, and fi is ith external force
at the location ri (i = 0, · · · , ∞). The potential ψ(x, t) is
given byψ(r, t) =∞
i=0

g 2D (r − ri ; t)fi ,(15)

where g 2D is the Green’s function in time domain.
The Fourier component of the potential ψ(r, t) is written by�(r, ω) =∞
i=0

G 2D (r − ri ; ω)fi ,(16)

where G 2D is the Green’s function in frequency domain.
In this paper, we define the Fourier transform as follows(Dahlen and Tromp 1998), ∞X(ω) =
x(t)e−iωt dt.
(17)∞

2.2 SI in an open system

Next, we consider SI in an open system. The formulationsin an open system depend on the source-receiver configuration: (1) random point sources distributed over thewhole space, (2) random sources distributed on a closed
curve, and (3) uncorrelated plane wave incidents fromvarious directions. This subsection describes the relationship between the different configurations of an open
system. Mathematically, the theory of an open systemdiffers from that of a closed system. We cannot use normal mode theory for an open system because the systemloses energy from the radiation boundary.
For ANT, CCFs are usually formulated in an open system, which better approximates the source-receiver configuration. Because the surface wave can be formulated
as a 2-D problem with a membrane approximation (e.g.,
Tanimoto 1990; Tromp and Dahlen 1992b, 1992a, 1993),
we can consider them as a 2-D potential problem in an
open system. For ANT, we explicitly express the mixedcomponent CCFs of surface waves in a simplified case at
the end of the subsection.
2.2.1 Seismic excitation by an infinite number of sources

We consider one realization of the background seismic
wave field excited by an infinite number of sources in an

We note that the Fourier convention depends on the literature. For example, Aki and Richards ’s definition has
the opposite sign in the exponential term.
Here, we consider a potential ψ(ro , t0 ) in a simplified
case with constant phase velocity as C(r) = C0, where rois the location of the origin and time denotes an arbitrary
positive time. The ψ(ro , t0 ) is represented by the sum of
the arrivals excited by sources along the concentric circle with radius r = C0 t0. The typical separation of thesources is assumed to be x . Within the circle with theband x , about 2π r/�x sources are distributed (Fig. 2left). Because the amplitude decay is proportional tor 1/2, the mean square amplitude ψ(ro , t0 ) is estimatedto be about 2π r/(�x(r 1/2 )2 ) = 2π/�x , which does notdepend on the distance r. Therefore, after t > 0, the fluctuations of ψ last for a semi-infinite time with the samemean squared amplitudes (Fig. 2 right).
2.2.2 Random sources distributed on a closed curve

The second source configuration is random sources on anarbitrary curve enclosing stations. Now, we observe thepotential ψ at ro within the circle with radius r (Fig. 2).
Based on the representation theorem, the wave excited by
distributed sources outside the circle can be completely
reproduced from the stresses and displacements on thecircle (e.g., Aki and Richards 1980). The contribution of

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

(a)Page 7 of 42

(b)

r
ro

t

∆x
t0

Fig. 2 Schematic figure for random sources in 2-dimension. a Source distribution in space. The red triangle shows the observed station. b Potentialψ(ro , t) observed at the origin r at time t. The wave packet at time t0 is traveled from the source within the concentric circle with radius r shownin (a)

the sources within the circle for the CCF can be neglected
because the contribution from the sources outside thecircle becomes infinitely large when considering infinitely
long times. These features mean that it is equivalent topersistent sources distributed only on the circle. Supposethat uncorrelated random sources with white spectra on
an arbitrary surface enclose stations in a heterogeneous medium. In this case, the time derivative of the corresponding CCF represents the exact Green’s function
between a pair of stations (Wapenaar et al. 2010a).
2.2.3 Uncorrelated plane wave incidents

The third source configuration is uncorrelated plane wave
incidents from various directions. Assuming that thesource-station distances are sufficiently longer than theaperture of the seismic array, the assumption leads to aplane-wave approximation (Nakahara 2006; Haney et al.
2012). The above discussions for a homogeneous source
distribution lead to the identical formulation of CCFs
regardless of the source configurations, as shown in thefollowing subsection. When we consider a heterogeneous
source distribution, the first source configuration is themost flexible to represent the source, and the second andthird are gradually less flexible. The flexibility of the firsttwo source configurations causes a complex dependenceof a CCF on the locations of the two stations. In contrast,
the inflexibility of the third source configuration causes a
simple dependence of the CCF only on the relative location between the stations (only the distance and the azimuth). In many cases, the ocean swell activities are far
enough away from the stations to approximate the phenomenon well. For simplicity, we will use the third source

configuration to represent CCFs based on plane wave
incidents in the following subsections.
2.2.4 CCFs for homogeneous source distribution in a 2‑D
homogeneous medium

For ANT, we explicitly express the CCFs between all
pairs of three-component seismometers in an open system when multimode Love and Rayleigh waves dominate
the ambient seismic wave field. We evaluated the CCFs
for the surface waves in a 2-D problem with a membrane
approximation (e.g., Tanimoto 1990; Tromp and Dahlen
1992a, b, 1993). We show an expression of the mixedcomponent CCFs of surface waves.
An arbitrary seismic wave field u(r, θ ; ω) in 2-D can be
represented by a superposition of multimode Love andRayleigh waves asu(r, θ ; ω) =

∞ 
∞


n=0 m=−∞

Ray

fn,m (ω)[Un (ω)Pm (r, θ ; ω)+Vn (ω)Bm (r, θ ; ω)]Love+ fn,m(ω)Wn (ω)Cm (r, θ ; ω),

(18)
Love
fn,m
where
is forcing of Love waves. Un and Vn are eigenfunctionsof the nth overtone of the Rayleigh wave, and Wn is the
eigenfunction of the nth overtone of the Love wave,
which has a real value. The basis functions Pm, Bm and
Wm are given byωreimθ ,
Pm (r, θ ; ω) = ẑJm Ray(19)
cn (ω)
Ray
fn,m (ω) is forcing for Rayleigh waves, andNishida et al. Progress in Earth and Planetary Science(2024) 11:4


 

Ray
ωr
θ̂ ∂
cn (ω) ∂r̂ +Jm Rayeimθ ,Bm (r, θ ; ω) =
ωr
∂r
r ∂θ
cn (ω)(20)
cLove (ω) r̂ ∂
ωr
∂
− θ̂
Cm (r, θ; ω) = n
Jm Loveeimθ ,ωrr ∂θ
∂rcn (ω)(21)
Ray
where cn (ω) is the phase velocity of nth overtone ofRayleigh wave, and cnLove (ω) is nth overtone of Love wave.
An assumption of equipartition of energy leads toRay∗ RayRay

�fn,m fn′ ,m′ �Un (ω)Un′ (ω) ≡ Pn (ω)δnn′ δmm′ ,(22)Page 8 of 42

Transverse
Radialr

1st station

2nd station

θ

RadialTransverse
Fig. 3 Schematic of the geometry of a pair of stations
with separation distance r and the coordinate. The figure shows
the radial and transverse directions for the station pair. To simplify
the equation, the signs of the directions are different from thoseof other studies. The difference changes the signs of CCF
in the components of RR, TT, RZ, and ZR

Love∗ Love
�fn,mfn′ ,m′ �Wn′ (ω)Wn (ω) ≡ PnLove (ω)δnn′ δmm′ , (23)
Ray

Ray

Ray∗

Love∗
fn′ ,m′ � = 0,
�fn,m fnLove′ ,m′ � = �fn,m(24)

Rn (ω) ≡ Vn (ω)/Un (ω),

(25)

Ray
where Pn (ω) is the power spectrum of nth overtone ofRayleigh wave, and PnLove is the power spectrum of nth

overtone of Love wave.
The mixed-component cross-spectra for multimode
Rayleigh waves are written as,

difference of sign RZ originates from the definitionof the radial direction, which differs from this study.
Using this formulation, cross-spectra for Rayleighwaves lead to a symmetric relation between the RZ andRayRay
ZR components as RZ = ZR , whereas cross-spectra
for body waves lead to an antisymmetric relation as
body
body
RZ = −ZR (Takagi et al. 2014).
Attenuation becomes significant when considering a seismic wave field in sediment above 0.1 Hz(Nishida et al. 2008a; Prieto et al. 2009, 2011). However,



�

�
�
�


|ω|r
ωr


−Rn J1 Ray
0 Ray Ray Ray  J0 cnRay

c


n




�
�
�
�
�ZZ �ZR �ZT
∞


�


2
Ray
 Ray Ray Ray 
R
|ω|r
ωr
n


,
0
P
(ω)
=
J
J
−R
 �RZ �RR �RT 
n
0−2
n 1
Ray
Ray


2
cn
cn




Ray
Ray
Ray
� 
�
n=0


�TZ �TR �TT





Rn2
ωr


0
0


Ray
2 J0+2(26)cnwhere a cross-spectrum represents the CCF in frequency
domain. The mixed-component cross-spectra for multimode Love waves are written as,physically plausible attenuation measurements arepractically difficult (e.g., Liu and Ben-Zion 2013; Tsai



0
0
0


Love
Love
�
�

∞ 
�Love


�
ZZ �ZR �ZT


1
ωr
Love
0
0 2 J0+2 cLove ,�Love �Love �Love  =
P
(ω)
n
n
RZ
RR
RT
�
�




Love
Love

1
ωr
n=0 
�Love


0
0
TZ �TR �TT
2 J0−2 cnLove

where J0−2 (z) ≡ J0 (z) − J2 (z) and J0+2 (z) ≡ J0 (z) + J2 (z).
Z represents the vertical component, and R and T represent the horizontal components according to the polarization direction (Fig. 3).
These equations are identical to the result of Haney
et al. (2012). A similar formulation for DAS observation is given by Nakahara et al. (2021). We note that the(27)2011) because there is ambiguity between the source
heterogeneities and the attenuation. Although newtechniques have been developed to overcome this problem (e.g., Liu et al. 2021; Magrini and Boschi 2021;
Bowden et al. 2015), they are beyond the scope of
this study. Even if the physically plausible attenuation
estimation is difficult, apparent Q measurements areNishida et al. Progress in Earth and Planetary Science(2024) 11:4

feasible for better fitting a synthetic cross-spectrum toan observation (e.g., Nishida et al. 2008a).
2.3 Travel‑time anomalies due to the source
heterogeneities

In this subsection, we evaluate travel-time anomalies due
to source heterogeneities. Although the mode decomposition in the previous subsection can be applied to moregeneral source distributions, the resultant CCF for heterogeneous source distribution becomes complex. Witha stronger assumption of the 2-D ambient seismic wave
field described by a superposition of uncorrelated planewaves from all azimuths (the third source-receiver configuration), the resultant CCF depends only on the relative location between the station pair. The condition thatthe excitation sources exist at infinite distances simplifies the problem (see Sect. 2.2.3 for details). Reducing
variables simplifies the effects on the source heterogeneities. This assumption is sufficiently realistic, because we
consider the on-land observation of the ambient seismic
wave field excited by distant ocean swells. This estimate isconsistent with Weaver et al. (2009), but with a different
derivation.
We assumed that the sources were located along aclosed curve surrounding the stations. If the source isdistant enough compared with the wavelength, the crossspectrum  between a pair of stations is given by�(r, θ ; ω) =∞im Jm (kr)βm (θ),

m=0where am and bm are Fourier coefficients of the source
distribution (Harmon et al. 2010; Cox 1973). The spatiallysymmetric part of the cross-spectrum with even m is areal function, whereas the antisymmetric part with odd
m is an imaginary function. The symmetric feature originated from the plane wave decomposition by the Bessel
function as
∞


m=−∞

Jm (kr) ≈

Nm (kr) ≈





mπ 
mπ  4m2 − 1
2−,sin χ −cos χ −
πkr
2
8kr
2

(31)



 


mπ 
mπ  4m2 − 12+,cos χ −sin χ −
π kr
2
8kr
2

(32)
where Nm is the Bessel function of the second kind (alsoknown as the Neumann function) and χ ≡ kr − π/4 . The
real part of the cross-spectrum can be evaluated by only
even orders as,∞Re[�(r, θ ; ω)] =

m=0

∞


Re[�(r, θ ; ω)] =(30)

Here, we evaluate a travel-time anomaly due to source
heterogeneities of CCFs by the asymptotic expansionof the Bessel function for a large argument given by
Eqs. 9.2.5 and 9.2.6 of Abramowitz et al. (1988),(33)(−1)m J2m (kr)β2m (θ)

m=0

2 sin χ ′′B ,πkr 2kr ev(34)
and the imaginary part of the cross-spectrum is given by≈ J0 (kr)Bev (θ) +Im[�(r, θ ; ω)] = i∞(−1)m iJ2m+1 (kr)β2m+1 (θ)

m=1≈ N0 (kr)iBod (θ) − i
where Bod and Bev are defined by

Bod ≡

Bev ≡



2 cos(χ) ′′B ,πkr 2kr od
(35)

∞


β2m+1 (θ),

(36)

∞


β2m (θ).

(37)

m=0

m=0

The real and imaginary parts lead to the following crossspectrum  as,(1)

�(r, θ ; ω) ≈ H0 (kr)B(θ)(2)

im Jm (kr)e−imθ .

(−1)m J2m (kr)β2m (θ).

Insertion of the asymptotic Bessel function into theabove equation leads to the following equations. The realpart of the cross-spectrum is given by(28)

where wavenumber k(ω) is ω/c(ω), m is the azimuthal
order, θ is azimuth (Fig. 4), and the Fourier coefficients ofthe source intensities βm (θ) are real value asa0
m = 0,2βm (θ) =am cos(mθ) + bm sin(mθ) m �= 0, (29)eikr cos θ =Page 9 of 42+ H0 (kr)B(θ + π )ieiχ B′′ (θ) − e−iχ B′′ (θ + π ) ,−
2π(kr)3(38)
where B ≡ Bev (θ) + Bod (θ). B(θ) represents the intensityof incident plane waves as a function of azimuth θ . Here,
we use the following relations.

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

y

B(θ)
δθ
r

θ

x

Fig. 4 Schematic figure of the geometry of stations
and the stationary phase regions shown in blue. The sources arelocated sufficiently distant from the stations, and the distribution
is given by B(θ ) as a function of the incident azimuth θPage 10 of 42

the higher frequency, the travel-time anomaly decreases
for a longer distance and a higher frequency. The equation does not include the first derivative of the intensity B
because the symmetry cancels the contribution (Weaver
et al. 2009). Equation 42 gives an error proportional to
the second derivative of the intensity B to the azimuth.
This equation can estimate the phase velocity bias due
to source heterogeneities and correction (Weaver et al.
2009; Froment et al. 2010).
Finally, the analytic representation of a cross-spectrum
for heterogeneous sources is given by π�(r, θ ; ω) ≈ 4 e− 2 i G ∗ (r; ω)e−iωδt(r,θ ) B(θ)
 (43)π+e 2 i G(r; ω)eiωδt(r,θ +π ) B(θ + π ) ,where Green’s function in 2-D is given by(2)

B(θ) + B(θ + π ) = 2Bev (θ),B(θ) − B(θ + π ) = 2Bod (θ).

G(r; ω) = −i

To evaluate the travel-time anomalies due to source
heterogeneity, we consider the far-field approximation(kr ≫ 1) of the Hankel function of the first kind:
12 i χ − 8kr
(1)
(39)
e
.
H0 (kr) ≈π kr
We also assume that the source heterogeneity is weak as

1 d2B
≪ 1.
krB dθ 2(40)

The causal part of the cross-spectra (Eq. 38) can be written byieiχ B′′ (θ )(1)
H0 (kr)B(θ ) − 
2π(kr)3
2
1
B′′ (θ )≈ B(θ )
exp i χ −
−
.
π kr
8kr
2krB(θ )(41)
Therefore, the travel-time anomaly δt of the causal partcan be estimated by
δt(r, θ) ≈B′′ (θ),2ωkrB(θ)(42)

where t is the travel time and ω is the nominal angularfrequency. Figure 4 shows the stationary phase regions,which dominate the contribution of the CCF (Snieder
2006). The aperture δθ of the stationary phase regions is
proportional to the square root of the ratio between the
wavelength and the station separation distance r (Froment et al. 2010). Because the narrower the aperture at

H0 (kr).
4(44)This subsection shows how to estimate travel-timeanomalies due to source heterogeneities. We note that
the CCFs still satisfy the original wave equation even if
the source distribution is heterogeneous. This property
ensures that phase velocity measurements can be possible even for a heterogeneous source distribution if thestation density is sufficiently high and the excitation
sources are located outside the station array. Following Lin et al. (2013), we briefly show the relation in thefollowing.
Here we consider potential ψ, which satisfies a wave
equation L as

 2
∂
2
− C(r)∇ ψ(r, t) = 0.
L[ψ(r, t)] ≡
(45)
∂t 2
We also assume that no sources are distributed within
the station array. The CCF φ is defined by
1 T
φ(r1 , r2 , t) = lim
ψ(r1 , τ )ψ(r2 , t + τ )dτT →∞ T 0(46)
The CCF satisfies the wave equation as
1
T →∞ T

L1 [φ(r1 , r2 , t)] = lim T0

L1 [ψ(r1 , τ )]ψ(r2 , t + τ )dt = 0,(47)where L1 is the wave equation concerning r1. The CCF
also satisfies the wave equation for r2. Because the CCF
satisfies the wave equation only without spectral whitening and one-bit normalization, these procedures often
create unphysical phases (Nakata et al. 2013).

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

3 Data processing
If the seismic wave field is stochastic stationary, dataprocessing for the CCF calculation is straightforward.
In a realistic situation, the spectral structure of seismic
records changes with time. Moreover, seismometers
record transient phenomena such as earthquakes, which
are noise for the CCF calculations. Instrumental noise isalso problematic. To improve the quality of CCF, we needto select and regularize seismic records (Bensen et al.
2007; Ritzwoller and Feng 2019). This section explains
the practical data processing.
3.1 Instrument

Before processing data, we must consider the influenceof instrumental characteristics. This review considers
the geophysical observable of ground motions and pressure in the ocean. The seismometer records the groundmotions, and they can be categorized into strong motion
sensors, short-period sensors, and broad-band sensors
according to the frequency and amplitude ranges.
3.1.1 Sensors

Today, on a global scale, broad-band seismometers arecommonly used. Because instrumental responses dependon sensor types, the responses should be removed prior
to data analysis. Typically, broad-band sensors have a flat
response in particle velocity above the corner frequency
of about 10 to 100 s, and the response decreases withdecreasing frequency below the corner frequency.
Short-period sensors (typically velocity seismometers with a corner frequency of 1 Hz) are also commonly used on a local scale. In the recent developmentof dense observation with more than 1000 stations (alsoknown as a large-N array, e.g., Lin 2013), geophones
with a natural frequency of 10 Hz are also often used.
Although only nominal responses are available in manycases, the responses are slightly different from one toanother (Ueno et al. 2015; Takeo et al. 2022). When we
use the data with a frequency close to the natural frequency, the difference might cause an apparent phaseshift, which could be problematic for a short separation
distance.
When a sensor is deployed on the seafloor, pressuregauges are also commonly used (Webb 1998). Cross-correlation analysis of ocean bottom pressure gauges exhibits
Rayleigh wave propagations by cross-correlation analysis(e.g., Takeo et al. 2014). When a sensor is deployed on theseafloor or bottom of a borehole, the orientation of thesensor is unknown (e.g., Takagi et al. 2019). In such cases,
the orientations should be determined by known events,
such as an active shot or teleseismic events.

Page 11 of 42

3.1.2 Instrumental corrections

In ambient noise cross-correlation analysis, small-amplitude seismic wave signals are used compared to regularearthquake signals. Moreover, the signals may fall within
the lower frequency band than the natural frequency ofthe seismometers, potentially residing below the instrument noise level. In such cases, instrument noise can significantly impact the quality of CCFs. Instrument noise
can be categorized into incoherent and coherent noises
across observation stations. Incoherent noise between
stations decreases with the square root of the number ofstacks of CCFs. However, coherent noise does not diminish and remains in the CCFs even with many stacks.
Takagi et al. (2015) conducted an ambient noise crosscorrelation analysis using Hi-net short-period seismometers in Japan and found periodic pulses at 60-secondand 1-second intervals within the CCFs. The stackingof raw records revealed that these periodic pulses were
not a result of data analysis but were inherent in the raw
records shared across all stations. Although the amplitudes were extremely small, less than 1 bit, they appearedin the CCFs at frequencies lower than the natural frequency band because of their coherence across stations.
Similar instrument noise and its impact on CCFs havealso been reported in other observation systems (e.g.,
Wang et al. 2017). This coherent periodic instrument
noise is attributed to load fluctuations in data loggers,
such as time calibration and data recording. There is also
temporally random but coherent instrument noise acrossstations. The coherent random instrument noise appears
as pulses near the zero lag in the CCFs. Such coherent
random noise was found in DAS observations (e.g., Tribaldos and Ajo-Franklin 2021) and cabled ocean-bottom
seafloor seismometers (Takagi et al. 2021).
It is crucial to remove coherent instrument noise toutilize CCFs to image subsurface structure. Several methods have been proposed to deal with periodic coherent
instrument noise. One method involves creating waveforms of periodic coherent instrument noise by stacking
many raw records and subtracting them from the raw
observation records (Takagi et al. 2015). Alternatively,
instead of computing instrument noise waveforms, one
can calculate the difference of raw observation records(e.g., the difference between two adjacent days) and thencompute the CCFs of these differential waveforms (Takagi et al. 2021). To reduce randomly coherent instrument
noise, practical methods such as modeling and estimating the spectral shape of instrument noise (Takagi et al.
2021) or subtracting the median of CCFs from each trace(Tribaldos and Ajo-Franklin 2021) have been proposed.
By taking advantage of dense DAS observations, the F-K

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

filtering of fast-propagating phases before computing the
CCFs is effective (Fukushima et al. 2022). In summary,
while coherent periodic or coherent random instrument
noise across observation stations can contaminate CCFs,
understanding the characteristics of the noise and applying appropriate removal methods enable us to expand theanalyzable frequency range and ultimately improve theresolution of structural imaging.
The time reliability is crucial for the cross-correlation
analysis between stations. For onshore observations,satellite systems, such as the Global Navigation Satellite
System (GNSS), usually give a precise time stamp. When
GNSS reception fails, the temporal change in the CCFs
has been used to estimate the error in the time stamps(Sens-Schönfelder 2008; Hirose and Ueda 2023). For offshore observations, precise time stamps are mostly given
at deployment and recovery, and the drift of the internal
clock can be linearly corrected. When the time stamp
cannot be obtained either at the beginning or at the endof the observation, the temporal change in the CCFs canbe used again to estimate the drift (Hannemann et al.
2014; Gouédard et al. 2014; Takeo et al. 2014). The timeasymmetric shape of the CCFs can also be used to estimate unknown instrumental responses such as the constant time shift of the logger or the frequency-dependent
phase response of the differential pressure gauges (Takeo
et al. 2014).
3.2 Data selection

We must choose seismic data that satisfy a local stationary state to apply the SI. Outliers such as earthquake
records and instrumental glitches decrease the SNR.
This subsection introduces (1) one-bit normalization, (2)RMS-based data selection, and (3) polarization-based
data selection.
3.2.1 One‑bit normalization

To suppress the effects of transient phenomena suchas earthquakes, one-bit normalization (e.g., Aki 1955;
Bensen et al. 2007; Cupillard et al. 2011) keeps only thesign of the original information, changing all positive
values to 1 and all negative values to -1. Because thismethod is simple, it has been widely used. When transient phenomena (e.g., many aftershocks and packet loss
during data logging) occur frequently, their contributions
decrease the quality of the CCFs. In such a case, careful reduction of transient signals allows us to improve
the CCF’s quality. A disadvantage of this method is theloss of amplitude information. Different amplitude

Page 12 of 42

normalizations among components cause problems in
the cross-correlation analysis of multicomponent data.
3.2.2 RMS‑based data selection

If the signal obeys a stationary process, we can reject outliers using an RMS threshold. In a strict sense, secondary
microseisms do not follow a stationary process but are ina local stationary state: RMS changes significantly by several orders of magnitude on a timescale of several days
(Fig. 5). We must define the local background levels forthe rejection, which change slowly over several days.
Here, we introduce an example for estimating the local
background level by Nishida and Takagi (2022). To define
the typical background level of the whole network at the
ith time step, we calculate the median of MS amplitudes
for all stations Pi . Here, we consider a situation: ith time
step is the latest accepted time step, and we reject n successive time steps. If Pi+n changes suddenly, we reject the(i + n + 1)th time step with the threshold ǫ:| ln(Pi+1 ) − ln(Pi )| > nǫ,(48)

we reject all data at time step i + n + 1 with the criterion
proportional to the rejection duration (orange points in
Fig. 5).
3.2.3 Data rejection associated with large earthquake using
a catalog

Seismometers record global propagations of many seismic phases excited by large earthquakes, typically with
moment magnitudes greater than 6 (Ekström et al. 2012).
In most frequency ranges, the RMS criterion can reject
the corresponding data. However, the secondary microseisms are still large enough to hide some earthquakes.
Even smaller amplitudes can bias the CCFs because theearthquake signals are coherent among the stations.
Careful rejections of large earthquakes improve the quality. Using a global earthquake catalog (e.g., the global
Centroid Moment Tensor (CMT) catalog (Ekström et al.
2012)), data rejection based on the magnitude in the catalog is feasible to exclude hidden seismic phases.
Approximately the amplitude of earthquake data
decays exponentially with time. The typical duration Deis:De = 2Q

ln(JM/1023 ),(2π f )

(49)

where f is the dominant frequency, Q is the typical quality factor, M is the moment, and J is a typical geometricalNishida et al. Progress in Earth and Planetary Science(2024) 11:4

Page 13 of 42

Fig. 5 Example of the data selection process. The orange dots show data excluded using the global CMT catalog (Ekström et al. 2012), and the blue
dots show data transients when the amplitude changes suddenly. The black dots represent selected segments. The vertical axis shows the relative
power normalized by the Peterson NLNM (Peterson 1993)), which represents the lowest ground noise on Earth. Taken from Nishida and Takagi(2022)

spreading (Nishida and Kobayashi 1999; Tanimoto andUm 1999). Data were rejected for De seconds after thearrival. The orange dots in Fig. 5 show a typical example
of rejected data (Nishida and Takagi 2022). Data can berejected until the earthquake signal is much smaller thanthe background noise levels.
3.2.4 Polarization‑based data selection

If the excitation processes are stochastic stationary intime and space, the energy partition among the mode
branches should be constant. The energy partitionchanges over time because they are not stationary ina realistic situation. For example, Takagi et al. (2018)
pointed out that the energy partition of P-wave microseisms at periods of 4–8 s becomes more significant on
seismically quiet days based on the polarization analysis of Hi-net data in Japan. When we want to emphasize overtones of seismic surface waves, data selection
based on polarization information can potentially
improve the detection of overtone branches. Pedersen
et al. (2023) proposed the data selection based on the
H/V spectral ratio to extract teleseismic body wave
microseisms, and this strategy may be applicable fordetecting overtone branches.
3.3 Regularization: weight and normalization of CCFs

If both the signal and the noise are subjected to theGaussian distribution, calculating the CCFs fromselected data is straightforward. However, neither noise

nor signal behaves well in a realistic situation. Regularization of data is required as part of the preprocessing(Bensen et al. 2007).
3.3.1 Weighting of cross‑spectra

The weighting on the data is important when we
calculate a cross-spectrum between a station pair
(Nishida 2014). The amplitudes of microseisms at
frequencies around 0.1 Hz change with time, whichreaches more than one order of magnitude on a timescale of one day. Spectral whitening efficiently reduces
the non-stationarity (e.g., Bensen et al. 2007). The
amplitudes of the seismic hum in the mHz band arestationary, although the local noise level is higher than
the signal levels. In this case, the weighting of the data
depending on the local noise level can be effective(Nishida 2014; Takeo et al. 2013). In the mHz band,the noise levels of the horizontal components areorders of magnitude higher than those of the verticalcomponents. For the calculation of the cross-spectrum, we suppressed noisy Fourier components using
the data weighting as follows.
We calculated a weighted cross-spectrum �ij (f )between the ith and jth stations at a frequency f as

�ij (f ) = k1
wijk (f )kwijk (f )Uik∗ (f )Ujk (f ),(50)

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

where Uik (f ) is a Fourier spectrum of ground accelerationof the kth segment at the ith station. The weighting factor wijk (f ) depends on the situation. If the noise is muchlarger than the signal, the weight can be estimated by

wijk (f ) =1

1|Uik (f )| |Ujk (f )|

.

(51)

We can evaluate the uncertainty σij (f ) of the resultant
�ij (f ) as

−1
Nij
�
1
1

σij (f ) = �
wijk (f ) ,(52)
Nij (f ) Nij (f )k=1

Page 14 of 42

where Nij (f ) is number of stacked traces (Takeo et al.
2013). The corresponding CCF φij (t) is calculated by the
inverse Fourier transform of �ij (f ) (e.g., Fig. 6).
3.3.2 Spectral whitening

If the signal is larger than the noise, the assumption of
the previous Sect. 3.3.1 is broken down. Although simple stacking works for signals with Gaussian distribution, the signal level changes significantly with time. In
particular, the amplitudes of secondary microseisms
can change by several orders of magnitude over a timescale of a few days. Without regularizing the amplitude,
the resultant CCFs emphasize days of high-ocean-swell
activities too much. For the amplitude regularization, theFig.
6 Cross-correlation functions as a function of interstation distances filtered at 20–30 s, 10–20 s, and 5–10 s. The amplitudes are multipliedwijk (f ) in the frequency domain to produce similar noise levels for different pairs of stations. The dashed pink line shows the group velocitybykof each specified mode at a typical period for the 1-D isotropic model. The modes include the fundamental mode and first-higher-mode Rayleighwaves, 0 S, and 1 S. The TT components show wave packets of Love waves xT, which could not be identified as a single mode due to possible
interference between the fundamental mode and the higher mode with very similar group and phase velocities. Taken from Takeo et al. (2016)

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Page 15 of 42

cross-spectrum ij is normalized by the amplitude (e.g.,
Aki 1957; Nishida et al. 2008a; Prieto et al. 2009) as,

�ij (f ) =kwijk (f )Uik∗ (f )Ujk (f ) = Uik∗ (f )Ujk (f )k|Uik (f )||Ujk (f )|.

(53)The ij with spectral whitening is also known as crosscoherency.The equation is similar to the weighting, but
the factor k wijk is not required.
Although direct waves are accurately reconstructedwith spectral whitening, they cause many pseudo-arrivals. Therefore, if we ignore the pseudo-arrivals, all normalizations work practically (Nakata 2020). To mitigate
the pseudo-arrivals, smoothing of wijk in the frequency
domain was proposed (Tauzin et al. 2019).
Another implementation is performed by spectral
whitening of the individual seismic trace before calculating CCFs (Bensen et al. 2007). Specifically, the inverse
Fourier transform of Uik (f )/|Uik∗ (f )| corresponds to thewhitened seismic trace.
3.3.3 Temporal flattening

Different amplitude normalizations at different stations
become problematic if we utilize amplitude information(e.g., attenuation). Although this topic is beyond the scope.
Temporal flattening (Zhou et al. 2020; Weaver 2011) was
proposed to preserve the amplitude information.

4 Dispersion measurement of surface waves
for constructing local 1‑D structures
In the first step of ANT (Fig. 1a), a reference dispersion
curve of phase or group velocities is measured using the
CCFs of an array. This section briefly summarizes themeasurement methods: (1) frequency–time analysis, (2)slant stack technique, (3) SPAC method, (4) FJ method,and (5) waveform fitting. The dispersion curves are also
important for the construction of local 1-D structures,
which can be an initial model of ANT (see Sect. 6.3.1).
4.1 Frequency–time analysis (FTAN)Frequency–time analysis is a common method for measuring the dispersion curves of earthquake data (Levshin
et al. 1992; Cotte and Laske 2002; Romanowicz 2020).
This method is also feasible for ambient noise CCFs
(Bensen et al. 2007; Harmon et al. 2007; Yao et al. 2011;
Spica et al. 2018). For the measurements, we usually use
the time-symmetric part of CCF to reduce the effects
of source heterogeneities (see equations 34 and 35 in
Sect. 2.3). In this subsection, we use the folded CCF tolag time 0 φij+ as,φij+ (t) = 12 φij (t) + φij (−t) t ≥ 0 ,
0
t<0

(54)

and the corresponding Fourier component is definedby �+ij (f ). The folding procedure distorts the waveformwhen the separation distance is shorter than about one
wavelength. The tail of the acausal wave packet contaminates the causal peak because the CCF does not satisfy
causality, as in the case of Green’s function.
A group velocity-period diagram can be calculated foreach CCF. The general outline of the procedures is as follows, although there are some other implementations.
Figure 7 shows an example of a synthetic test.
At a given center frequency f, a narrow-band filter withthe center frequency f0 is applied to the CCF. For example FTAN (Levshin et al. 2018; Ritzwoller and Feng 2019)use a Gaussian filter

Fb (f , f0 ) = e

−α



f −f0
f0

2

,

(55)

where α is the coefficient, determining the bandwidth(Fig. 7a, b). Because a smaller α broadens the bandwidth,a smaller α provides a stable estimate for noisy data. At
the same time, smaller α can bias the estimation in casesof strong dispersion. α often changes with the separation
distance (Ritzwoller et al. 2011) as
α = α0 r/r0 ,(56)

where α0 is the reference value at the separation distance
r0. In a synthetic test shown in Fig. 7, α0 is 40, and r0 is
200 km.
The envelope functions of the filtered CCFs are plotted against the periods and corresponding group velocities calculated from the separation distance and the lagtime (Fig. 7c). Correction for f0 using the instantaneous
frequency based on the analytical signal is feasible for abetter estimate (Shapiro and Singh 1999; Levshin et al.
2018). In this case, a synthetic CCF is given by�syn (ω) = J0 (k(ω)r) + n(ω),(57)

where n is Gaussian noise with a standard deviation of
0.03. We can measure group velocities by tracing the lociof a local maximum as a function of the period. Because
a strong frequency dependence on the spectral content
can bias the group velocity measurements, the correction is necessary for accurate measurements. AFTAN(Automatic Frequency–Time Analysis) package (Barmine2018) is available for this type of measurement. We show
an example of the estimation errors at 0.1 Hz as a function of the separation distance (Fig. 7d).

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Page 16 of 42(a) Gaussian filter

1.0

10km
20km
50km
100km
200km
500km

0.5

0.0
0.00

0.05

(b) Gaussian filter

1.0

0.10

0.15

0.5

0.0

0.20

0.0

0.1

Frequency [Hz](c) FTAN diagram

30.5(d) Anomaly of phase velocity and group velocity

2.8Group velocity
Phase velocity

2.6
2.4

1
0
−1

2.2
2

4

6

8

10

12

14

16

18

−2

20

0

200

Period [s]∆= 3 λ180(e) Phase

Number density

0

− 90

200

400

600

800

1000

(f) Errors of Group/phase velocity

90

0

400

Distance [km]5Phase [degree]

0.4

2

3.0

− 180

0.3

∆= 3 λ

3.2

Anomaly [%]Group velocity [km/s]

3.4

2.0

0.2

Frequency [Hz]

600

800

1000

Distance [km]Phase velocity

4
3
2Group velocity

1
0
− 1.0

− 0.5

0.0

0.5

1.0

1.5

2.0

Anomaly [%]Fig. 7 a Gaussian filters for various separation distances at 0.1 Hz. b Gaussian filter for r = 200 km. c FTAN diagram for a synthetic CCF
with the observation error is 1%. The separation distance is r = 200 km. The black line shows the theoretical group velocity and the red dots show
measured values using synthetic data. d Group and phase velocity anomalies for the synthetic CCFs with 1% measurement errors as a functionof separation distance at 0.1 Hz. e Measured phases (blue dots) and the theoretical values (dashed lines) at 0.1 Hz. f Histogram of group and phasevelocities estimation errors at 0.1 Hz

We can also measure the phase velocity from thephase of �+ij (f ) (e.g., Levshin 2018; Ritzwoller and Feng2019) if a single mode is dominant. Multimode interference is problematic for such measurements. To isolate a single mode, we must apply a proper filter. Here,

we introduce two typical filters: (1) time variable filter(Landisman et al. 1969) and (2) floating filter (Levshin
et al. 1992).
The time variable filter (Harmon et al. 2007; Landisman et al. 1969) isolates a single-mode wave packet by

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Page 17 of 42

windowing the CCF around the arrival time predictedby the group velocity measurements. Window length isset according to the dispersion because the stronger dispersion makes the duration of the wave packet longer. A
typical window length (Harmon et al. 2007) Tw is given by

Tw = 160 ×

dU
,
dτ

(58)

where U is group velocity, and τ is period.
Floating filtering is also a common procedure forextracting a single mode (Levshin et al. 1992). Based on
the group velocity measurement, phase correction ψ(ω)
is defined by ω1
ψ(ω) = rU −1 dω + t1 ω,(59)
ω0where r is the separation distance, and t1 is a constant.
iψ�+ij (ω) is multiplied by e to compress the signal aroundthe average group time. t1 is chosen to focus on the average group arrival. After cleaning up the signal, φ + isrecovered by the inverse Fourier transform with e−iψ.
The phase velocity C can be measured from the phaseof �+ (ω) (Lin et al. 2008; Harmon et al. 2007) as,C(ω) =

ωr− arg(�+ (ω)) + 2π N + π4

,

(60)

where N is an integer number. The π/4 originates fromthe Bessel function, and the 2π ambiguity can be solved if
a good starting model is available. Note that �+ (ω) is outof phase by π/2 with the corresponding empirical Green’s
function (see Sect. 2 of Lin et al. (2008)).
Figure 7e shows the phase of the synthetic test and thetheoretical values. Although they agree, the finite bandwidth of the Gaussian filter biases the measured phases.
At higher frequencies, estimating the phase ambiguity N
is difficult due to the ambiguity of an initial model. If adense array is available, the array methods, as discussed
in the following subsections, are feasible for measuring the phase velocity at higher frequencies. The measurement using FTAN has a limitation for the separation
distance of CCFs. The group velocity measurements areaccurate when the separation distance exceeds three
wavelengths. For this reason, (Bensen et al. 2007) recommended that the longest period should be r/12.
Histograms in Fig. 7f show the estimation errors of
the group and phase velocities. Although the measured
phase velocities were relatively accurate, the dispersion
caused bias. When applying this method to strongly dispersive wave fields, we must carefully consider bias, particularly for low-SNR data.
An advantage of the FTAN method is the ability tomeasure the dispersion of a single station pair. As alreadyFig. 8 Group and phase velocities and the average spectrogramof raw data. Gray-scale filled contours represent the averaged
spectrograms of all the CCF that met the SNR requirements,
where the maximum value for each period in the individual
spectrograms was normalized to one. The fundamental and second
mode Rayleigh wave is visible in these data as the broad regions.
Filled black circles with error bars indicate the average groupvelocities for the fundamental mode, while open circles with error
bars indicate the average group velocities for the higher-mode
Rayleigh wave. Black triangles with error bars indicate fundamental
mode uncorrected phase velocities; gray triangles with error bars
indicate −π/4 corrected phase velocities for higher and fundamental
mode phase velocities, and open triangles with error bars indicate
higher mode phase velocities. Taken from Harmon et al. (2007)

pointed out, the FTAN method also has the disadvantagethat the dispersion measurements are difficult for separation distances shorter than about three wavelengths.
FTAN method is still feasible for measuring multimodedispersion if a wave packet of each mode branch is isolated in the time domain. If multiple phases overlap intime domain, an array-based method (slant stack or
SPAC) is appropriate.
Figure 8 shows an example of a multimode Rayleighwave dispersion diagram measured using ocean-bottom
seismometers. The figure shows two clear branches. The
fundamental mode with a period 8 s has energy in theocean, while it has energy in the crust and the uppermost mantle with the longer period. The overtone also
has energy in the crust and the uppermost mantle. Thisbehavior of a multimode Rayleigh wave in an oceanic
region is common, as shown in Fig. 6. A seismic arrayon a soft sediment layer records multiple mode branches
above about 0.2 Hz (e.g., Spica et al. 2018).
4.2 Slant stack technique

Suppose that a dense array with station spacing comparable to the wavelength is available. In this case, the slant
stack technique, also known as the frequency–wavenumber (F-K) method, is a common array processing (RostNishida et al. Progress in Earth and Planetary Science(2024) 11:4

and Thomas 2002; Gouédard et al. 2008). This method
also assumes that the superposition of plane waves canrepresent the wave field. The beam power Bp is definedby a sum of time-shifted waveforms �+ (ω) with expectedtravel-time anomaly ωp · xi

2
 √


+
iωp·(xk −xl ) Bp (p, ω) = 
rkl �kl (ω)e
(61)
 ,


kl

where k and l represent the station number, rkl is theseparation distance between kth station and lth station,and xk is the location vector of kth station represented bya Cartesian coordinate with the origin at the array cen√troid. The term rkl corrects the amplitude decay due togeometrical spreading for a surface wave.
The disadvantage of this method is a potential bias dueto the plane wave assumption. Equation 31 shows that
the J0 (vertical component) the J0−2 /2 (radial and transverse component) can be approximated as√π1πrJ0 (kr) ≈ cos(kr − ) +sin(kr − ),
4
8kr
4

(62)

√
π
π
r
7
J0−2 (kr) ≈ cos(kr − ) −sin(kr − ). (63)
2
4
8kr
4
The contribution of the terms (8kr)−1 and 7(8kr)−1 is significant for a shorter station pair. These equations suggest that the bias of the horizontal components could
be larger. We tested the bias of the vertical component
using a synthetic CCF at 0.15 Hz given by Eq. 57 and theslowness of 0.3 s/km. Figure 9 shows the station distribution of the array (Hi-net stations with a radius of about
100 km) and the result. The beam peaks at 0.3 s/km, consistent with the theoretical value. To verify the bias, wemeasured the peak slowness for 10,000 experiments. Figure 10 shows the histogram. Although the measurements
fluctuate, the central value is 0.05% lower than the theoretical value. The bias of about 0.05% is consistent with−(8k 2 r 2 )−1 estimated from equation 62 with a typical
separation distance of 60 km.
This method can reveal multimode dispersion on ascale of 100 m to 100 km: DAS observation on a 100 m
scale (Dou et al. 2017), a 10 km scale ocean bottom DAS(Viens et al. 2022; Williams et al. 2021), and a 100 km
scale basin structure (e.g., Boué et al. 2016; Jiang and
Denolle 2022). Li et al. (2020) demonstrated that supervised machine learning methods are feasible for separating multimode information.
There is another implementation to calculate the F-K
spectrum. When calculating the slant stack, the separation distances between stations are not regular. We canPage 18 of 42

re-sample the CCFs on regular grids if all the station separation distances are shorter than half the wavelength. CCFs
can be mapped from the spatial domain to the F-K domainusing the 2-D fast Fourier transform (FFT) after a correction r 1/2 for geometrical spreading (Gabriels et al. 1987).
This method is applied to joint inversion of the first overtone and fundamental mode for deep imaging in the Valhall oil field using ambient seismic noise (Tomar et al. 20172018).
4.3 SPAC method

This subsection explains the implementation of the SPAC
method. In the case of homogeneous source distribution
synin a stratified medium, the synthetic cross-spectrum ρZZ
between vertical components can be represented by a Bessel function (see Sect. 2.2.4) as ωr syn,ρZZ (r, a, C; ω) = a(ω)J0(64)
C
where r is the separation distance of the station pair.
Assuming an arbitrary wavenumber k(ω) at a given frequency, we determine the optimum amplitude a, by minimizing the squared difference S between the observed
cross-spectrum i of ith pair and ρ syn (ri , a, C; ω) (e.g.,
Nishida et al. 2008b) with weight wi . S is given by

S(a, C; ω) =

N

i=1

 syn
2
wi ρZZ (ri , a, C; ω) − �i (ω) , (65)where N is number of the pairs. By minimizing S at agiven frequency ω , we infer the phase velocity C(ω) andthe amplitude a(ω). Regarding a, S can be minimized
analytically as

∂S
= 0.
∂a

(66)

We obtain the optimum a for given C and ω asNwi �i (ω)J0

aopt (C, ω) = i=1N

i=1w i J0 ωr 
i

C

(67)

.

 ωr 2i

C

To minimize S with respect to C in ω , we maximize thevariance reduction VR given byN

VR (C, ω) = 1 − i=12 ωr i
− �i (ω)
wi aopt (ω)J0
C
N

i=1

.

2

wi �i (ω)(68)

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Page 19 of 42(a) Station distribution(b) Histogram of separation distance

500

38°N

400

Counts

37°N

36°N

300

200

35°N
100

34°N
135°E

136°E

137°E

138°E

139°E

0

140°E

0

50

100

Distance [km]

150

200

(c) Plot in slowness domain
1.0

SPAC
Slant stack
FJ

Amplitude

0.5

0.0

0.1

0.2

0.3

0.4

Slowness [s/km]
Fig. 9 Synthetic test of slant stack, SPAC, and FJ methods. We assumed an ambient seismic wave field with a slowness of 0.3 s/km. a Station
distribution. b Histogram of separation distance between station pairs. c Plot in slowness domain: slant stack, variance reduction of the SPAC
method, and the FJ spectrum as a function of slowness. The noise is assumed to be 3% of the signal

SPAC: VR

Number density

20

Beamforming
FJ: Hu et al. (2023)
100

FJ: Nimiya et al. (2023)

-0.1

0

0.1

0.2

Slowness anomaly [%]Fig. 10 Bootstrap errors of the estimated slowness by SPAC,
slant stacking, and FJ method. The noise is 3% of the signal,and the sampling number for the bootstrap method is 10,000

0.4

0.5

The maximum is inferred using the grid search method.
The estimate can be refined using a generalized leastsquares method (Menke and Jin 2015).
Figure 9 shows a test of the SPAC method using a
synthetic cross-spectrum at 0.15 Hz given by Eq. 57
and the slowness of 0.3 s/km. We used an array with aradius of about 100 km shown in Fig. 9a. Most of theseparation distance ranges from 40 to 130 km (Fig. 9b).
Figure 9c shows VR as a function of slowness (1/C).
Because the amplitude of the surface wave decreases√
with 1/ r , the appropriate modeling by SPAC sharpens
the peak compared to the beam power.
We measured the maximum of VR for 10,000 experiments to evaluate the estimation errors. Figure 10
shows the histogram. The measurements fluctuate
within 0.05%, and the median value is around the theoretical value.

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Modeling RR and TT in horizontal components ismore complicated because both Love and Rayleigh waves
contribute to these components, as shown by Eqs. 26 andsyn
27. The corresponding synthetic cross-spectra ρRR andsynρTT can be represented by� ωr �� �� R �
� syn �
ωr
J
J
0−2
0+2
ρRR
R�
cL � a (ω)
c

�
�
.
=
syn
(69)
ωr
aL (ω)
ρTT
J0+2 R J0−2 ωrLcc

When we analyze the station with a separation distanceshorter than the wavelength, we must consider the Love/
Rayleigh wave for the RR/TT component, respectively.
Although the Love and Rayleigh waves are coupled inthis system, we can solve this equation directly. The
maximum search becomes complex because the variancereduction at a given frequency is a function of C R and C L.
Here, we consider a simplified problem with a farfield approximation. Because J0+2 is proportional tor −3/2, the term becomes negligible for r longer than the
wavelength (e.g., Takeo et al. 2013). In this case, we cansimply apply the SPAC method to these components
separately by replacing J0 of the vertical component with
J0−2. Figure 11 shows VR using the SPAC method in Japan(Nishida et al. 2008a). The diagram shows clear Rayleigh and Love branches. Although this method assumed
a single-mode branch, the figure exhibits multimodebranches, which enable us to increase depth resolution(e.g., Ikeda et al. 2012). The SPAC method is also applied
to seismic data from the ocean floor (Takeo et al. 2013,2014; Lin et al. 2016; Takeo et al. 2016, 2018; Kawano
et al. 2023). To separate multimode Rayleigh waves,polarization information could also be useful (Nayak andThurber 2020).

Page 20 of 42

On a global scale, cross-spectra can be modeled using
the Legendre function instead of the Bessel function(Nishida et al. 2002; Nishida 2014).
4.4 FJ method

The frequency–Bessel (FJ) method (Wang et al. 2019a;
Hu et al. 2020; Nimiya et al. 2023) is closely related to
Aki’s SPAC method. Here, we reinterpret the FJ methodwith the SPAC method. For simplicity, we consider the
ZZ component in this subsection. Refer to Hu et al.
(2020) for other components.
If we observed φ(r, ω) as a function of the separation
distance r, the squared difference with weight r can beevaluated by the following integration,
 rN 
 ωr 2
− �(r, ω) rdr.
S=a(ω)J0(70)
C
0
By minimizing S for given C and ω , we obtain rN  
ωr
J0�(r, ω)rdr
C
aopt (C, ω) = 0 rN   2
ωr
rdr
J0
C
0

(71)

If rN is much longer than the wavelength, the denominator can be approximated by rN   
ωr 2
rN C
J0
.
rdr ≈
(72)
C
πω
0

aopt (C, ω) can be approximated by the integral (thenumerator) evaluated at discrete points (Hu et al. 2020)
asFig. 11 Frequency–phase velocity spectra using all pairs of stations. a A spectrum of radial components shows a clear Rayleigh wave branch,and b that of transverse components, which shows a fundamental Love wave branch, first overtone branch, and second overtone branch. Taken
from Nishida et al. (2008a)

Nishida et al. Progress in Earth and Planetary Science

aopt (C, ω) ≈(2024) 11:4

 ωr πω i�i (ω)J0r̄i �rirN CCi ωr  r 2 + 2ri (ri+1 − ri−1 ) − r 2
πω 
i
i−1
i+1
�(ri )J0
,
=
rN C
C
8
i

(73)

where



rj + rj+1
1 rj−1 + rj
+
,
r̄i ≡
2
2
2
�r ≡

1
(rj+1 − rj−1 ),
2

(74)

(75)

and i is the cross-spectrum of the ith station pair. FJ
spectrum I(C, ω) defined by Nimiya et al. (2023) can beapproximated by ωr 
rN ω opt
ω2 rN
�(r, ω)J0
rdr ≈a (C, ω).
I(ω, C) ≡
C 0
C
π
(76)
Although Hu et al. (2020)’s definition of the coefficientω2 /C 3 instead of ω2 /C has a physical meaning, theemphasis on small phase velocity C may bias the estimation, as shown below. We can interpret the FJ spectrum
as a variant of SPAC method with weight wi defined bywi =2
2 + 2r (rπ ω ri+1
i i+1 − ri−1 ) − ri−1
.
rN C
8

(77)

We again tested the FJ method using a synthetic crossspectrum given by Eq. 57 and the slowness of 0.3 s/km.
Figure 9c shows the FJ spectrum. The peak width is narrower, but the side lobes are more significant than theothers. Let us evaluate the side lobes using Eq. 71. If rN is
sufficiently longer than the wavelength, the FJ spectrum
ω for synthetic data J (ωri /C0 ) can be approximated (seeEq. 31) by


 ωr
ωri
π
π
2ω 
i
−
−
�ri .
cos
cos
I(C, ω) ≈
π
C
4
C0
4
i
(78)
This integral can be interpreted as the Fourier series.
Fourier series without a taper function cause significant
spectral leakage (Oppenheim and Schafer 2014). If rN is
sufficiently large, I(C, ω) has value only at C = C0 likethe delta function, due to the orthogonality of the cosine
function (more precisely, the orthogonality of the Bessel
function). This means that the FJ method has, in principle, the highest resolution in phase velocity.
To evaluate the accuracy, we measured the peak slowness for 10,000 experiments using the FJ method. Figure 10 shows the histogram. The estimated value ofthe FJ spectrum of Hu et al. (2020) is about 0.2% larger

Page 21 of 42

than the theoretical value. The bias is caused by the falloff of the spectra at low phase velocity of C −2 (Shapiro and Singh 1999). The drop in the C −2 term (Nimiya
et al. 2023) of equation 78 significantly reduces bias, as
shown in Fig. 10. The greater variation in the estimatesis due to the greater weight of the smaller amplitudes atlonger separation distances (Eq. 77). When the observation error is sufficiently small, the high resolution of the
FJ method in the slowness domain will be beneficial for
separating multimode surface waves. When the observation error is significant, the SPAC method can consider
the noise using appropriate weights.
Figure 12 shows an example of field data (Hu et al.
2020). Although this subsection mentioned only the ZZ
components, it can be applied to other components.
Because the RR and TT components record both Love
and Rayleigh waves, as shown by Eq. 69, their coupling
complicates the estimation of the FJ spectra in the RR
and TT components. Hu et al. (2020) proposed procedures to separate the Love and Rayleigh wave using
orthogonal relations of the Bessel functions (Fig. 12).
4.5 Multimode dispersion measurements by waveform
fitting in a model space

The frequency–phase velocity diagram does not show
clear mode branches in regions with strong lateral heterogeneities, such as volcanoes and ocean sediments(Takeo et al. 2022; Takagi and Nishida 2022). For example,
Fig. 13 shows an example of SPAC diagrams in an oceanicarea offshore the central part of northern Japan near theJapan Trench. They show fundamental Rayleigh and Lovewaves below 0.08 Hz, but other mode branches are intermittent and ambiguous. Lateral variations of the shallow
soft sediment distort the dispersion curves of multiple
modes. Moreover, multimode interference also complexes
the diagrams. In this case, the physically plausible phasevelocity measurements are difficult. A physically plausible
constraint on the dispersion is crucial for better measurements. One strategy is to skip phase velocity measurements using a SPAC method. For a given S-wave velocity
structure, the corresponding phase velocity of jth overtone C j (β; ω) can be evaluated, where β is a vector withcomponents of S-wave velocity at each depth. The variance reduction given by Eq. 68 is summed over in a frequency range from ω0 to ω1 and the N mode branches as

VRall (β) =N−1  ω1
j=0

ω0

VR (C j (β; ω); ω)dω,(79)

where VRall is the summed variance reduction. By maximizing VRall for β, we can infer the local 1-D structure. This method can be applied to the RR and TT

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Page 22 of 42

Fig. 12 Multicomponent frequency–Bessel spectrograms of field data in North America. a–d are the FJ spectrograms IR0, IR1, IR2, and IL0, respectively.
Here, IR0 is FJ spectrum of Rayleigh wave using ZZ components, IR1 is FJ spectrum of Rayleigh wave using RZ and ZR components, IR2 is FJ spectrumof Rayleigh wave using RR and TT components, and IL0 is FJ spectrum of Love wave using RR and TT components. The black arrow in (a) indicatesthe location where the fundamental mode Rayleigh wave is bifurcated into two branches. The ith overtone is labeled next to each dispersion curve,and 0 indicates the fundamental mode. Taken from Hu et al. (2020)

Fig. 13 Array-based dispersion measurement of the real data. a 1-D shear velocity structure modeled with six layers of variable thickness(parameterization A). The thick and thin black lines represent the mean and standard deviation of the estimated models from the 100
bootstrap samples. b Rayleigh wave dispersion curves estimated with model A. The thick and thin white curves represent the mean and threetimes the standard deviation of the 100 bootstrap results. The background image is the average of the variance reductions of the conventionalcross-spectral fitting for the ZZ, ZR, and RR components. The frequency range without the white shadow is used in the dispersion measurement. c
Love wave dispersion curves estimated with model A. Taken from Takagi and Nishida (2022)

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

components even if Love and Rayleigh waves appear inboth the RR and TT components with shorter separation distances. Figure 13 shows an example. Figure 13ashows the S-wave velocity models, and Fig. 13b, c shows
the corresponding mode branches. Along the calculated mode branches, VR is summed over the frequency
and all wave types. By maximizing VRall , we can infer the
local 1-D model. This method searches for the global
maximum of VRall in the model space. Physically implausible models can be rejected automatically. The plot of VR
against frequencies and phase velocities, as shown in the
figures, is feasible to check the validity of the estimatedmodel. Liu et al. (2023) proposed an inversion method,which directly compared multimode dispersion diagramswith the synthetic kernel of the Green’s function, whichshares a similar idea as the method here.

5 Dispersion measurements for each path
In the second step of ANT (Fig. 1b), phase or groupvelocity is measured for each path. This section introduces four methods to measure dispersion for each CCF:
(1) FTAN, (2) the zero-crossing method, (3) the SPACbased waveform fitting in the phase velocity domain, and(4) the waveform fitting in a model space. FTAN and thezero-crossing method (Ekström et al. 2009) are feasible tomeasure the phase velocity of a single mode. The latter
two methods measure phase velocities by SPAC-based
waveform fitting. The first measures the phase velocity
directory under an assumption of single-mode dominance, whereas the second measures by fitting synthetic
CCF to the observed one with the help of reference CCF
constructed by the local 1-D structural inversion.
5.1 FTANFTAN explained in Sect. 4.1 can also be used for phaseand group velocity measurements of a single trace (e.g.,
Lin et al. 2008). We can measure the group velocities ofmultimode surface waves by picking up local maxima
from a stacked FTAN diagram over an area (Spica et al.
2018). However, multimode interference makes phasevelocity measurements using FTAN difficult. Even in thecase of multimodes, if we can isolate a single mode by atime variable filter (Landisman et al. 1969) or (2) a floating filter (Levshin et al. 1992), we can measure the phasevelocities.
5.2 Zero‑crossing method

The zero-crossing method (Ekström et al. 2009) is briefly
summarized here. We will consider the ZZ component
for simplicity, but it can be easily extended to other components. For the frequency ωn of the nth zero-crossing
observed, the phase velocity can be evaluated byPage 23 of 42

Cm (ωn ) =

ωn r
,
zn+2m

(80)

where zn of is the nth zero of the Bessel function J0,
m = 0, ±1, ±2, which represents the ambiguity of n.
This formulation allows missed or extra zero crossings
because observed noise makes measurement of zn difficult. The messed or extra zero crossings 2m arise fromthe distinction between positive-to-negative zero andnegative-to-positive zero. The advantages of this method
over the FTAN method are (i) applicability even for ashorter separation distance (typically shorter than three
wavelengths) and (ii) smaller bias for shorter distances
owing to no far-field approximation. Although zerocrossing measurements are practically robust (Ekström
et al. 2009), the low-SNR CCF behaviors become unstable. In particular, because strong lateral heterogeneities
distort the phase information, appropriate filtering prior
to application is crucial for stable measurements.
5.3 SPAC‑based waveform fitting in phase velocity domain

This method measures the phase velocity anomaly
in the ith frequency band from ωi to ωi+1 under theassumption of single-mode dominance. If we alreadyinferred local 1-D models based on the SPAC method,
the synthetic cross-spectrum can be represented byωr,�synth (ǫ; ω) = a(ω)J0(81)
C(1 + ǫ)where ǫ is phase velocity anomaly. We infer δa and ǫ byminimizing the squared difference defined by ωi+1 2S(αi , ǫi ) ≡(1 + αi )�synth (ǫi ; ω) − �obs (ω) dω,ωi(82)
where ǫi is the phase velocity anomaly at ωi , and αi is anamplitude correction term at ωi . αi is estimated first by
the least-squares method for ǫi . With the estimated αi , ǫiis measured using the grid search method (Nagaoka et al.
2012; Yamaya et al. 2021).
5.4 Multimode dispersion measurements by waveform
fitting in a model space

In the case of the dense spacing of the mode branches
in the phase velocity domain, there is a method for estimating the multimode dispersion curve by waveform
fitting using a 1-D velocity structure as model parameters. Yoshizawa and Kennett (2002) developed theoriginal idea of measuring the phase velocity dispersion curve for each path of teleseismic surface waves.
We applied this method to the cross-spectra of ambient

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Page 24 of 42

noise continuously (Takagi and Nishida 2022; Takeo et al.
2022). The advantages of this method are that it provides
physically achievable dispersion curves and can simultaneously measure dispersion curves for multimodesurface waves. Given the 1-D velocity structure, it is possible to calculate dispersion curves for all corresponding higher-order modes, explaining observed waveforms
through the superposition of multimode surface waves.
For the superposition of multiple modes, it is necessary
to know the excitation amplitudes for each mode. In thecase of teleseismic surface waves, the excitation amplitudes for each surface wave mode can be calculated fromthe seismic source mechanism. We can use the amplitudes aj calculated from the array-based SPAC method
for the cross-spectrum of ambient noise (Eq. 67). Given areference 1-D structure of the first step (Sect. 4), the synthetic cross-spectrum can be represented by

�

synth

(δβ; ω) =

N
−1

j=0

j

a (ω)J0ωr, (83)C j (β + δβ; ω)where δβ is a perturbation of 1-D structure, N is the number of the mode branches, and j represents the order ofthe mode branch. We infer δβ by minimizing the residualsum of squares given by2 ω1 w(ω)2 �synth (δβ; ω) − �obs (ω)S(δβ) =
dω.
w(ω)2
ω0
(84)
It is worth noting that the 1-D velocity structure, which
serves as a model parameter, does not necessarily represent the true Earth structure (Yoshizawa and Kennett 2002). The 1-D structure model is not unique withrespect to band-limited dispersion curves; different 1-D
models can reproduce similar dispersion curves. The
main products here are the multimode dispersion curves
calculated by the 1-D structure. Thus, the search range ofthe model space should be wide enough to reproduce a
dispersion curve that nicely fits the waveform.
These equations correspond to the vertical component but can be readily extended to other components.
In waveform fitting, a relatively wide frequency range isoften handled, and multicomponent data can be simultaneously processed. Therefore, it is important to consider the weighting of the data in practical applications.
The error estimate of the SPAC method based on thearray can be utilized as the weights (Takagi and Nishida2022).
Figure 14 exemplifies its application to S-net data
(Takagi and Nishida 2022). Even in this case, where
the signal-to-noise ratio of the CCFs is low, it is possible to robustly estimate multimode dispersion curves

Fig. 14 Dispersion measurement for a single path by waveform
fitting in a model space. a Modeled pair-specific 1-D shear
velocity structure for the S-net station pair N.S2N08-N.S3N20.
The thick and thin black lines represent the mean and standard
deviation of the 1-D structures estimated by the bootstrap
method. The gray line results from the reference model obtainedby the array-based SPAC method in Fig. 13a. b Phase velocity
dispersion curves of the fundamental-mode Rayleigh wave (0R),
the first-overtone Rayleigh wave (1R), and the fundamental-mode
Love wave (0 L). The thin blue and red lines are the means and threetimes standard deviations of the estimated pair-specific dispersion
curves. The thick blue and red curves represent the frequency
ranges that meet the quality control criteria based on variancereductions within narrow frequency bands. The gray ones are
the results of the array-based measurement in Fig. 13b, c. c Fittingof the cross-spectra for this pair. The black and red are the data
and the modeled cross-spectra, respectively. d Same as in (c),but both cross-spectra are normalized by the standard deviationof the cross-spectral data estimated from the array-based SPAC
method. Taken from Takagi and Nishida (2022)

by utilizing the entire waveform of multicomponent
CCFs and applying physical constraints based on 1-D
structure.

6 2‑D multimode ANTIn the third step of ANT (Fig. 1c), we conduct a 2-D
inversion of the phase/group velocity of multimode surface waves from the phase/group velocity measurements
for all pairs of stations. The method has been developed
to explore the mantle structure using earthquake data
(e.g., Woodhouse and Dziewonski 1984; Nakanishi andAnderson 1982; Tanimoto and Anderson 1985; Ekström
et al. 1997; Nolet 2008). Because review papers covering
ANT are available (e.g., Barmin et al. 2001; Montagner,

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

2015 Rawlinson et al. 2010; Thurber and Ritsema 2015; ),see them for details.
For inversion, we should consider the sensitivity kernels of the phase/group velocity. Although in an idealized situation, the kernel is identical to a finite-frequency
kernel for an earthquake, a heterogeneous distribution of
the noise source distorts the kernels (e.g., Nishida 2011;
Tromp et al. 2010; Hanasoge 2013; Fichtner 2014; 2015;
Fichtner et al. 2016; de Vos et al. 2013).
The first subsection describes how source heterogeneities affect the kernel using an analytical method for a
simplified case. The next section explains the conditions
for the application of surface wave tomography. The last
subsection briefly summarizes the 2-D inversion methods: local 1-D inversion, ray-theoretical inversion, andfinite-frequency inversion.
6.1 2‑D phase velocity sensitivity kernel for ANT

The sensitivity kernel for 2-D phase velocity is evaluatedbased on the Born and Rytov approximations (Born et al.
1999; Ishimaru 1997; Nolet 2008). The Born approximation relates the scattered wave to the phase velocity perturbations with the first-order approximation, but it does
not separate the effects on the amplitude and the phase.
On the other hand, the Rytov approximation by taking
the logarithm of the wave (e.g., Yoshizawa and Kennett

K (r1 , r2 , r3 ; ω) ≈Page 25 of 422005) enables us to separate the logarithmic amplitude
term and the phase term.
Here, we consider the Born sensitivity kernel K, which
relates the phase velocity perturbation and the crossspectrum perturbation δ� as
δC
δ�(r1 , r2 ; ω) =K (r1 , r2 , r3 ; ω)d�3 ,(85)
Cwhere C is phase velocity and δC is the perturbation. The
locations r1, r2 , and r3 are shown in Fig. 15a. The kernel K
is given byK (r1 , r2 , r3 ; ω) = − 2k 2 �(r2 , r3 ; ω)G ∗ (r1 , r3 ; ω)+�∗ (r1 , r3 ; ω)G(r2 , r3 ; ω) ,(86)
where �(r1 , r2 ; ω) is the cross-spectrum between r1 andr2 at frequency ω (e.g., Nishida 2011), and G is the Green’s
function in frequency domain. Although the kernel was
generally calculated numerically, it can be estimated byobserved data in principle if dense data were available(Chmiel et al. 2018).
With an assumption of the 2-D ambient seismic wave
field described by a superposition of uncorrelated planewaves as in Sect. 2.3, Eq. 43 gives the analytic representation of . Assuming that ωδt and kr are small enough, thekernel can be rewritten by−ik  +i(χ23 +χ13 )B(θ2 ) − e−i(χ13 +χ23 ) B(θ1 )e√
π r13 r23
+ei(χ13 −χ23 ) (B(θ2 + π ) − B(θ1 + π )) ,

(a)

(87)

(b)
Station 1

Station 2

θ1

r 13 > r 23

θ2

r 12
r1

r2

r 13

r 13 < r 23

r1

r2

r 12
r3

Fig. 15 a Schematic map of the geometry of stations. The star symbols show the station locations at r1 and r2. The open circle shows the locationof a phase velocity anomaly at r3. b Elliptic curves with constant |r13 − r23 |, and hyperbolic curves with constant |r13 − r23 |Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

where B(θ) represents the intensity of incident planewaves as a function of azimuth θ , and rij is the distance
between ri and rj , χij ≡ krij − π/4 (see Fig. 15a). The
form of a phase sensitivity kernel for a homogeneous
source distribution is the same as that of the earthquake
data (Nishida 2011).
The Rytov approximation is feasible to evaluate thephase sensitivity kernel (e.g., Yoshizawa and Kennett 2005; Nishida 2011). Although the approximation
requires a single-wave packet, the cross-spectrum 
given by Eq. 43 (θ = 0, and r = r12) includes both thecausal and acausal parts. We divide the cross-spectrum
 and the kernel K into causal and acausal parts according to the sign of the exponent, approximately. The causal
part can be represented by e−kr+ωt , while the acausal partcan be represented by ekr+ωt , where the wavenumber k
and the angular frequency ω are positive. The cross-spectrum  is divided into the causal part c and the acausal
part ac as,(2)�c (ω) = H0 (kr12 )eiωδt(r12 ,π ) B(π )(1)�ac (ω) = H0 (kr12 )e−iωδt(r12 ,π ) B(0).

(88)

The Born sensitivity kernel K is also divided into thecausal part K c and the acausal part K ac according to thesign of the exponent as,

Most studies on ANT used the time-symmetric part ofthe cross-correlation function to measure phase velocity¯
anomalies. The time-symmetric part of cross-spectrum 
is defined by

¯
≡
�(ω)

�c + �ac∗
,
2

where

¯
δ �(ω)≡

δ�c + δ�ac∗
2

K ≡

ac



(95)

K a (r1 , r2 , r3 ; ω) + K ac∗ (r1 , r2 , r3 ; ω).
2
(96)
The Rytov approximation is employed to obtain aphase sensitivity kernel for phase velocity perturbations (Nishida 2011). In the Rytov method, the loga¯ is considered instead ofrithm of the cross-spectrum K̄ (r1 , r2 , r3 ; ω) ≡





(93)

¯ and
We consider the relation between perturbation of 
the corresponding Born kernel K̄ as

δc
¯ =
δ�
K̄ (r1 , r2 , r3 ; ω)d�3 ,(94)c−ik(χ13 +χ23 ) B(θ ) + eik(χ13 −χ23 ) (B(θ + π ) − B(θ + π ))
√−ik
1
2
1
π r13 r23 −e −ik(χ +χ )−k13
23 B(θ )√−e1
π r13 r23

cPage 26 of 42

r13 < r23
r13 ≥ r23

,

(89)

and

K

≡

ik(χ13 +χ23 ) B(θ )√−ik
r13 < r23
2
π r13 r23 e


,
−k
ik(χ
+χ
)
ik(χ
−χ
)
13
23
13
23
√
B(θ2 ) + e(B(θ2 + π ) − B(θ1 + π )) r13 ≥ r23
π r13 r23 e



where the region r13 > r23 corresponds to the first andfourth quadrants, whereas the region r13 < r23 corresponds to the second and third quadrants (Fig. 15b).
The Born kernels also satisfy K = K a + K ac from thedefinition.
The causal perturbation δ�a and acausal perturbationδ�ac are related to phase velocity anomalies as

δc c
c
δ� =
K (r1 , r2 , r3 ; ω)d�3 ,
(91)
c

δ�

ac

=



δc ac
K (r1 , r2 , r3 ; ω)d�3 .
c

(92)

(90)

¯ can
the cross-spectrum itself. By taking the logarithm, be divided into real (amplitude) and imaginary (phase)parts, as follows:¯ = ln(Ae−iψ ) = ln A − iψ,
ln �

(97)

where A is the amplitude of the causal part, and ψ is its
phase. The phase perturbation δψ can be written by

δc
δψ =
Kp (r1 , r2 , r3 ; ω)d�3 .
(98)cwhere the phase sensitivity kernel Kp is the imaginary¯.
part of −K̄ /�

Nishida et al. Progress in Earth and Planetary Science(a) The first term

(2024) 11:4

Page 27 of 42(b) The second term(c) The whole kernel

x3Distance [km]

100
50

r2

r1

0
− 50

− 100

− 100 − 50050Distance [km]100

− 100 − 50050Distance [km]−10Amplitude

100− 100 − 50050Distance [km]1001Fig. 16 a The first term of the phase kernel Kp for B = cos θ + 1. The first term has an elliptical shape. The yellow triangles show the station location.
b The second term of the phase kernel Kp. The second term has a hyperbolic shape. c The whole kernel of the phase kernel Kp. The backgroundphase velocity is 3 km/s. The frequency range is from 0.1 to 0.2 Hz. We also apply a cosine taper for averaging in the frequency range

Then, the phase kernel Kp can be written by

k 3/2
Kp = √
2π





r1
π  B(θ2 ) + B(θ1 )− cos k(r13 + r23 − r12 ) −r13 r23
4 B(0) + B(π )π  B(θ2 + π ) − B(θ1 + π ),+ sin k(|r13 − r23 | + r12 ) +4
B(0) + B(π )where we neglect the first-order term.
Figure 16 shows an example of sensitivity kernels. The
separation distance of the station pair is 150 km, and thephase velocity is 3 km/s. The source B(θ) is assumed tobe cos θ + 1.
The first term has an elliptical shape (Fig. 15b), as
shown in Fig. 16a, corresponding to the finite-frequency
kernel of an earthquake. The kernel along the line
between the station pair (θ1 = 0 and θ2 = π ) is identical to the finite-frequency kernel of an earthquake. The
width of the kernel at the midpoint is proportional to√r12 , where  is the wavelength 2π/k . As a location
moves away from the line between, the source heterogeneity increases the sensitivity on the r2 side.
The second term with the dependency of |r13 − r23 |exhibits a hyperbolic shape (Fig. 15b) caused by source
heterogeneity (Fig. 16b), and the term oscillates rapidly
outside the stationary zones behind it. The effect tends tobe smoothed out when the typical spatial scale of the lateral heterogeneity is larger than the wavelength.
When the spatial scale of the source distribution islarge (i.e., the azimuthal dependence of B(θ) is smooth),the shape of the kernel is similar to the finite-frequency(99)

kernel of an earthquake (e.g., Spetzler et al. 2002; Yoshizawa and Kennett 2005). Because an approximation of
the finite-frequency kernel of an earthquake works formost cases practically, most studies did not consider thesource heterogeneity effect of the kernel. To check theeffects, travel-time anomalies due to the source heterogeneities given by Eq. 42 can be a common criterion. In
particular, because source heterogeneities cause apparent
azimuthal anisotropy, the travel-time anomalies given by
Eq. 42 are useful to evaluate accuracy. Although attenuation measurements are beyond the scope, the corresponding amplitude kernel is more sensitive to the source
distribution, so that the source heterogeneities should beconsidered.
6.2 Conditions for application of surface wave
tomography

This subsection briefly summarizes the conditions for
applying the 2-D inversion method, known as phase/group velocity tomography. Surface wave tomography requires weak scattering of the surface waves during propagation because we measure phase or group

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Page 28 of 42

Fig. 17 Log–log plot of reported isotropic scattering coefficients against frequency. The left vertical axis shows the isotropic scattering coefficientgiso, and the right vertical axis shows the corresponding mean free path ( g−1
iso ). Grey lines are for the lithosphere. Taken from Fig. 1 of Sato (2019)velocity anomalies of the direct waves. Because the propagation distance ranges from hundreds to several thousand km, the distance should be shorter than the mean

Fig. 18 ka − kL diagram (Aki and Richards 1980) for a classificationof wave propagation, where k is the wavenumber, a is the scaleof heterogeneities, and L is the propagation distance. Red lines show
propagation distances from 10 to 150 km, and blue lines show thosefrom 150 to 2000 km. The mean free path at 0.03 Hz is 5000 km,
that at 0.1 is 500 km, and that at 0.2 Hz is 30 km

free path. Otherwise, strong scattering makes phase orgroup velocity measurements difficult. Figure 17 shows
the mean free path as a function of the frequency compiled by Sato (2019). When we explore the mantle structure using earthquake data, we use surface waves below
0.05 Hz. The mean free path longer than 1000 km at thefrequency validates the application of the method of surface wave tomography.
To understand applicability, Fig. 18 shows a classification of wave propagation and the applicable method inka − kL diagram (Aki and Richards 1980), where k is thewavenumber, a is the scale of heterogeneities, and L isthe propagation distance. Here, we assume that the mean
free path is comparable to a, although this estimate isvery crude. The mean free path of northern Japan from
0.5–1Hz is estimated to be about 30 km (Hirose et al.
2020), and that of Germany is estimated at 0.2 Hz to beabout 500 km (Sens-Schönfelder and Wegler 2006b). The
red lines show propagation distances from 10 to 150 km,
and the blue broken lines show those from 150 to2000 km. At 0.03 Hz, the surface wave propagation at teleseismic distance beyond 2000 km is within the regime of
ray theory. At 0.2 Hz, the strong scattering makes measurements difficult for propagation distances longer than
200 km. This is the reason why earthquake surface wave
tomography is practically difficult above 0.1 Hz. Before

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4(a) Earthquake surface wave tomography

Page 29 of 42(b) ANT

Fig. 19 a A typical example of earthquake-receiver geometry for earthquakes. The gray lines show ray paths between the earthquakes
and the stations. Red stars show the earthquakes, and black dots show stations. b A typical example of station geometry for ANT. The lines show ray
paths between the station pairs

applying ANT, setting up an appropriate problem basedon the observed situation is important.
Let us consider the difference in ray geometry between
ANT and earthquake surface wave tomography. In thecase of regional-scale earthquake surface wave tomography with source-receiver paths, the target area is surrounded by earthquakes, which are limited along active
tectonic areas (Fig. 19a). Because most ray paths travel
across the area, the travel-time anomalies reflect theseismic structure along lines over the scale. Because theperturbations of the inferred phase/group velocity map
strongly depend on the damping of the tomographic
inversion, the absolute velocities tend to be ambiguous. Although the damping problem is common in raytheoretical inversion, including ANT, an initial local 1-D
model for ANT constructed from multimode dispersion
measurements using subarrays (see Sect. 6.3.1) can mitigate the ambiguity.
In contrast, in the case of ANT with a dense array,the distribution of the ray path is uniform (Fig. 19b).
Above 0.05 Hz, CCFs with longer separation distances(typically longer than 1000 km) become complex. Ray
paths with separation distances shorter than 1000 km
enable us to infer tomographic maps even in regionswith strong lateral heterogeneities above 0.1 Hz. Thissituation is similar to earthquake tomography using
the two-station method (e.g., Dziewonski and Hales
1972; Hamada and Yoshizawa 2015), which measures

the phase differences between seismograms for two stations on a common great-circle path. The drawbacksof the two-station method are (1) that it requires both
stations to be located close to the common great-circle
path, which reduces the available ray paths, and (2) thatthe longer earthquake-receiver distance tends to cause
the wave propagation to enter the regime of multiple
scattering at a higher frequency (Fig. 18).
The shape of the sensitivity kernel is elliptical along
the ray path,√ and the width at the midpoint is proportional to r12 , where  is the wavelength. Because thelengths of ray paths of earthquake data become longer
in general, the widths become wider than those of
ANT. The wider width tends to average the anomalies
within the fat ray, and the averaging over a long distance makes restoring the small-scale image difficult.
The advantages of ANT are the homogeneous path distribution and the shorter ray paths. These features enable
us to estimate phase/group velocity maps above 0.1 Hz.
6.3 Implementations of the 2‑D phase/group velocity
inversions

The subsection briefly summarizes implementations of
2-D phase/group velocity inversions using phase/groupvelocity anomalies measured from ambient noise CCFs:(1) regionalization using subarrays, (2) the ray-theoretical inversion of ANT, and (3) finite-frequency inversion
of ANT.

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Page 30 of 42

6.3.1 Regionalization using subarrays

If dense stations are distributed over a large area, the stations can be divided into local subarrays. Phase/groupvelocity maps can be inferred from phase/group velocity measurements with the subarrays, such as F-K analysis (Jiang and Denolle 2022), FJ method (Fu et al. 2022;
Nimiya et al. 2023), and slant stack (Qin et al. 2022).
These methods can easily incorporate information on
multimodes with smaller amplitudes because array analysis enhances the SN ratios.
The inferred map is also helpful as an initial model
for further inversions. Most ANT studies employ linearized inversion to obtain 2-D phase/group velocity maps (e.g., Rawlinson and Sambridge 2004) because
the computational cost of fully nonlinear inversion (e.g.,Markov Chain Monte Carlo (MCMC); Bodin and Sambridge 2009; Bodin et al. 2012; Galetti et al. 2016) is stillhigh. A good initial model is crucial for better modelingfor linearized inversion, because an inappropriate initialmodel may lead to a solution with a local minimum. An
inappropriate initial model sometimes requires strong
damping to stabilize the inversion. The appropriate initialmodel could mitigate the degradation of a tomographic
image.
6.3.2 Ray‑theoretical inversion of ANT

The most common method of tomographic inversion ofphase/group velocity anomalies is ray-theoretical inversion with the assumption of a high-frequency limit. For
example, tomo_sp_cu_s is a program (Barmin 2018)
designed for estimating 2-D models from regional- or
global-scale surface wave phase/group velocity measurements (Barmin et al. 2001; Ritzwoller 2002). The program
also assumed a straight ray for a weakly heterogeneous
medium. Another implementation to solve the eikonal
equation is the fast marching method (Sethian 1996;
Rawlinson and Sambridge 2004; White et al. 2020). SURFTOMO (Rawlinson 2008) is an iterative nonlinear traveltime tomography code in 2-D spherical shell coordinates.
Figure 20 shows an example of wavefront tracking using
the final tomographic results. Although the figure shows
the concentric wavefront within the arc distance of 10◦ , itshows kinks in Eastern Australia, which can be explainedby triplication due to the low-velocity zones. This figure
also suggests that CCFs with a separation distance longer
than 10◦ are difficult to model above 0.1 Hz without thehelp of a good initial model. An advantage of ANT is that
we can choose ray paths that meet the conditions according to the problem setting owing to the homogeneous
distribution of the ray paths.
The eikonal tomographic method (for example, Linet al. 2009) also assumes a high-frequency limit with asingle mode and solves the eikonal equation directly. The

Fig. 20 a Evolution of the wavefront in the case of strongly
heterogeneous media. The station with a green star acts
as the source, and the red triangles show the stations
concerning the advancing wavefronts. Note the triplication effects
around the low-velocity zones. b Histograms of passage time
residuals for measurements at 10 s for Rayleigh waves, initiallycompared to a uniform reference model with a velocity of 2.8 km/s
with great-circle paths, and after the inversion with 6 iterationsof the subspace method with deviated paths. Taken from Fig. 4of Saygin and Kennett (2012)

eikonal tomographic method tracks the evolution of the
wavefront using only local information with dense arraydata. The phase velocity can be estimated by the in situ
spatial gradient of the phase travel time without aninversion. The method can be applied to a more generalcondition of the source distribution than equal energy
distribution because the method does not require CCFs
to be related to the corresponding Green’s function. Thismethod works as long as the phase travel time satisfies
the eikonal equation (see the end of Sect. 2.3; Snieder
et al. 2006; Lin et al. 2013; Lehujeur and Chevrot 2020).
If a dense array satisfies the station spacing shorter than
the wavelength, this method is feasible to infer the phasevelocity map (e.g., de Ridder and Dellinger 2011; Lin et al.
2013; Ritzwoller and Feng 2019; Chen et al. 2022).

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

6.3.3 Finite‑frequency inversion of ANT

Although the patterns of the ray-theoretical inversions
are robust in different studies, the retrieved strength ofphase/group velocity perturbations differs (e.g., Yoshizawa and Kennett 2004). Finite-frequency inversion is
crucial to improve the amplitude recovery of velocity perturbations. As shown in Sect. 6.1, the sensitivity kernel
has an elliptical√
shape, and the width at the midpoint is
proportional to r12 . For a heterogeneous data set withmultiple scales, the difference in the width of the sensitivity kernels becomes problematic, because it becomes difficult to recover the amplitude of velocity perturbations
by a ray-theoretical inversion. This subsection summarizes recent developments on finite-frequency inversions.
As explained in Sect. 6.1, a finite-frequency sensitivity
kernel of a CCF also depends on the source distribution.
Consequently, the phase/group velocity structure has a
trade-off with the source distribution (Fichtner 2015). To
improve the structural image, we must infer the source
distribution simultaneously (Fichtner et al. 2016; Sager
et al. 2020, 2017). To mitigate the effects of source heterogeneities, measurements of differential travel times
between two pairs of three stations were proposed (Liu2020). The sensitivity kernels of the three-station method
significantly reduce the effect of the heterogeneous
source distribution. Because CCF amplitudes are moresensitive to source heterogeneities (Tsai 2011), the considerations of the source information will become morecrucial in inferring the attenuation structure.
Because the computational cost of the sensitivity kernels for both the structure and the source distribution
is still high, the earthquake sensitivity kernel for phase/group velocities is instead used practically (Gao andShen 2015a; Yang and Gao 2018; Covellone et al. 2015;
Gao and Shen 2015b; Emry et al. 2019; Yang et al. 2022;
Gao and Shen 2014; Wang et al. 2019b 2020; Lu et al.
2020; Crosbie et al. 2019; Janiszewski et al. 2019; Russell
and Gaherty 2021; Zhao et al. 2020). The approximation
becomes exact under the condition of equipartition ofenergy (Nishida 2011). Although the kernel does not consider the source heterogeneities, it improves the image
from ray-theoretical inversion.
Another implementation of the finite-frequency effectof ANT is Helmholtz tomography (Lin and Ritzwoller
2011b; Lin et al. 2012b; Mordret et al. 2013). This method
is a natural extension of eikonal tomography to consider
finite-frequency effects by introducing the amplitude
correction term. Although this method requires dense
data, the phase velocity can be estimated without thesensitivity kernel and the inversion, as in the case of eikonal tomography. This method can be applied to a moregeneral condition of the source distribution (see the endof Sect. 2.3).

Page 31 of 42

7 3‑D seismic velocity inversion using
the multimode dispersion
In the fourth step of ANT (Fig. 1d), a 3-D seismic velocity structure is constructed from the inferred phase/group velocity maps in the previous step, which exhibit
the dispersion curves at each grid point. The local 1-D
structure at each grid point can be inverted from thedispersion curves. 3-D seismic velocity structures canthen be constructed by collecting local 1-D structures.
In most cases, only the dispersion of fundamental modes(Love and Rayleigh waves) is utilized for the inversion
because the excitation sources are distributed on the surface. Recent studies reveal that a dense array enables us
to utilize multimode dispersion. Joint inversion with the
H/V or receiver function is feasible to improve the depth
resolution.
7.1 1‑D seismic velocity inversion at each grid

The measured dispersion curves are inverted to obtain
local 1-D structures. Inversion methods can be categorized into a Monte Carlo method (Shapiro and Ritzwoller
2002; Maraschini and Foti 2010) and a linearized method(Tarantola and Valette 1982; Tarantola 2005; Montagner 2015). Because the local 1-D inversion is identical tothe earthquake data, this section briefly summarizes therelated studies mentioned above.
All 1-D inversions require a code to calculate the surface wave dispersion (e.g., DISPER80, Saito 1988; senskernel-1.0, Levshin 2018; Computer Programsin Seismology, Herrmann 2013). These codes also
calculate the depth sensitivity kernel, which is also useful for a Monte Carlo inversion. Figure 21 shows a typicalexample of depth sensitivity kernels for a model based on
a reference Earth model PREM (Dziewonski and Anderson 1981), where the thickness of the ocean is changed
from 3.0 to 4.6 km, and the thickness of the crust is
reduced from 22 to 6 km (Takeo et al. 2013). Fundamental Love and Rayleigh waves have sensitivity at shallower
depths. The higher modes are not sensitive to the S-wave
structure at shallower depths, and they are more sensitive
to the S-wave velocity structure at greater depths thanthe corresponding fundamental models. With increasing periods, the kernels have higher sensitivity at greater
depths. A multimode inversion can improve depth resolution because overtones have different sensitivities fromfundamental modes.
Figure 22 shows an example of the 1-D inversion(Yamaya et al. 2021) using multimode dispersion under
the ocean. The study used a computer program package
SEIS_FILO (Akuhara 2022) for the inversion, which
implements a transdimensional Bayesian interface andreversible-jump MCMC (rj-MCMC) (Green 1995). They
assumed a uniform distribution for the prior probability

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Page 32 of 42

(a)

(b)
6

1st Lovec (km/s)

5
4
3

1st

h

leig

Ray

Depth (km)00th Love
0th Rayleigh

2
1

5

10

20

50

5

0

0th Rayleigh

10s
5

10

10

0.0

0.5

1.0

1st Rayleigh5s

0th Love

0

10s

5s

0.0

0.5

1.0

Kα (km−1)

0

5s

10s

Depth (km)15

Kα (km−1)0

1st Rayleigh10s

25s 5sPeriod (s)

(c)

0

5s

15

100

0th Rayleigh

1st Love

50

10s
100

10s
25s

150

25s

5s

25s
200
50

0.00

0.05

Kβ (km−1)

0.10

50

0.00

0.05

0.10

Kβ (km−1)

50

0.00

0.05

Kβ (km−1)

0.10

0.000

0.005

Kβ (km−1)

Fig. 21 a Phase velocities of the fundamental (0th) and first-higher modes of surface waves for the modified PREM (Dziewonski and Anderson
1981). The sudden phase velocity decrease marked by the ellipse corresponds to the transition of modal energy from solid to ocean. b, c
Sensitivities of phase velocities to P-wave and S-wave velocities (Kα and Kβ , respectively). The two lines in each panel show the depths of the seafloor
and Moho. Taken from Takeo et al. (2013)

distribution (1–10 for the number of layers; 2.3–10 km
for depths of the discontinuities; and 2.3 km for seafloor
depth and 0.1–5 km/s for S-wave velocities). P-wave
velocities and densities were scaled by the empirical
relation (Brocher 2005). The multimode dispersion constraints enable us to infer both sediment thickness andS-wave velocity. To demonstrate the feasibility of multimode measurements to improve depth resolution, let
us compare the inversions of Fig. 22 with and withouthigher-mode surface wave dispersion. Figure 23d shows
the inversion result using multimode dispersion, whereas
Fig. 23e shows the result using only fundamental modes.
These figures show that the inversion using multimodedispersion is feasible to constrain the shallow sedimentstructure, improving the resolution of the crustal structure. Although Rayleigh waves also have sensitivity to
the P-wave velocity, the limited sensitivity only at shallow depths (Fig. 21b) makes it difficult to constrain the
P-wave velocity structure practically. With an empirical
relation (e.g., Brocher 2005), the P-wave velocity and thedensity are scaled to the S-wave velocity to reduce thenumber of parameters.
Many groups have already published crust and uppermost mantle models using fundamental modes, although

this paper focuses on multimode ANT. Since there is nospace to discuss individual models, we introduce Data
Services Products: EMC, A repository of
Earth models (Trabant et al. 2012; IRIS DMC 2011),
which is a resource that offers access to different Earthmodels, along with visualization tools to preview models, facilities to extract data and metadata from the models, and access to software and scripts that others in the
research community have contributed.
Seismic anisotropy gives us a clue to understanding therheology and deformation of the Earth. In some cases,
an isotropic S-wave velocity model cannot explain theobserved discrepancy between the Love wave dispersion
and the Rayleigh wave dispersion, which can be explainedby radial anisotropy (e.g., Aki and Kaminuma 1963;
Anderson 1962). Radial anisotropy in the upper mantle
can be explained by lattice-preferred orientation in the
low-velocity zone or partial melt under shear stresses(e.g., Montagner 2015). Radial anisotropy also exists inthe lower crust, which can be explained by the latticepreferred orientation (LPO) of anisotropic crustal minerals under extension (Moschetti et al. 2010a, b; Huang
et al. 2010). ANT also revealed radial anisotropy in volcano areas (Jaxybulatov et al. 2014; Harmon and Rychert

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Page 33 of 42

Fig. 22 a Stations of the ocean bottom seismometers array in the region off Ibaraki. The station intervals are about 6 km. The yellow triangles show
the stations deployed on October 17, 2010, and the red triangles show those deployed on February 15, 2011. b, c S-wave velocities along the A-A’
and B-B’ lines. The red, orange, and green regions reflect layer 1, layer 2, and layer 3, respectively. We interpret the red and orange regions to be
sedimentary layers, the green region to be Cretaceous sediment, and the blue region to be the upper crust, respectively. The black lines show
the depth of the acoustic basement (between the orange and green regions) and the top of the upper crust (between the green and blue regions).
Modified from Figs. 1 and 8 of Yamaya et al. (2021)

2015; Spica et al. 2017; Lynner et al. 2018; Jiang and
Denolle 2022; Jiang et al. 2023, 2018; Miller et al. 2020).
The layering of a partial melt layer in the magma reservoir can explain them. Thus, radial anisotropy is feasiblefor discussing the stress state and the texture. Azimuthal
anisotropy can also provide information on rheology anddeformation (Lin et al. 2009; Yao et al. 2010; Ritzwoller
and Feng 2019; Russell and Gaherty 2021): for example,
the flow in the uppermost mantle related to plate motion(Takeo et al. 2016). We note that apparent anisotropy
due to the lateral heterogeneities of the isotropic S-wave
velocity structure is problematic (Lin and Ritzwoller
2011a; Fichtner et al. 2013). When estimating anisotropy
by the tomographic method, the anisotropy has tradeoffs with the lateral heterogeneities. Moreover, ANT has
another difficulty that source heterogeneity can yield up
to about 1% apparent azimuthal anisotropy (Harmon
et al. 2010; Takeo et al. 2014, 2016), which can be estimated by Eq. 42 (equivalent to Weaver et al. (2009)).

7.2 Joint inversion with other observations

Although multimode ANT improves depth resolution,
surface wave tomography has an inherent low depth resolution. Joint inversions with other geophysical observables are feasible for a better depth resolution at the 4thstep.
A thin low-velocity layer near the surface is problematic. Rayleigh waves are essentially less sensitive to theS-wave velocity structure near the surface, as shown in
Fig. 21. The poor sensitivity of Rayleigh waves at shallow depths makes resolving the shallow structure difficult, although the layer can still change the dispersion
curves of surface waves. The polarization of the Rayleighwave provides information on the shallow S-wave velocity structure beneath a seismic station. If the amplitude
ratios between the horizontal and vertical components(ellipticity) for individual mode branches are available,
the joint inversion can constrain shallower depths (Linet al. 2012a). Even if the ratios of individual modes are40.2
Frequency(b) First overtone

3
2
1
0.2
Frequency

Phase velocity (km/s)

0.1

4

0.3

0

(c)Fundamental Rayleigh wave
0.1

3
2
1
0.1

0

0

(e) Fundamental mode

0.1

layer 1layer 2layer 3

5

10

0

1
2
3
4
S wave velocity (km/s)

5

0

0.05

5

10

0

1
2
3
4
S wave velocity (km/s)

5Posterior probability

Phase velocity (km/s)

0.10(d) Multi-mode

Depth (km)1Page 34 of 42

Posterior probability

2(2024) 11:4

0

Depth (km)3

Posterior probability

0.1

Posterior probability

4 (a) Fundamental Rayleigh wave

Posterior probability

Phase velocity (km/s)

Nishida et al. Progress in Earth and Planetary Science00.2Frequency

Fig. 23 Comparison of 1-D S-wave velocity structure inversions using both the fundamental and the first-higher modes of Rayleigh wave (a,b, d) with that using only the fundamental modes (c, e). a Phase velocity of the fundamental mode of the Rayleigh wave. The blue points show
the average phase velocities. The Markov chain Monte Carlo method calculated the posterior probabilities for 1-D average S-wave velocity
structures. b Phase velocity of the first-higher mode of the Rayleigh wave. Blue points show the 1-D average phase velocities. c The phase velocityof the fundamental mode of Rayleigh waves with posterior probabilities using only fundamental modes. d S-wave velocity structure inferredby the MCMC method using both the fundamental and the first-higher modes of Rayleigh wave. The blue line shows the median velocity at each
0.1 km depth grid point. Layers 1 and 2 are well constrained despite the given loosely bounded uniform priors. e S-wave velocity structure inferredby the MCMC method using only the fundamental modes. Modified from Figs. 5 and 10 of Yamaya et al. (2021)

difficult to measure, the spectral ratio of the horizontal
to vertical (H/V) components at a station (Nakamura
1989; Sánchez-Sesma et al. 2011; Bahavar et al. 2020)
can be robustly measurable. The H/V ratio allows us toconstrain the shallow structure by joint inversion with anassumption of equipartition of modes (Spica et al. 2018).
We note that the ellipticity also has information on theazimuthal anisotropy (Lin and Schmandt 2014) because
the ellipticity also depends on the azimuth from one station to the other station.
In particular, single-mode surface wave inversion is difficult to distinguish between the smooth Vs model andthe layered model because the dispersion curves of surface waves are insensitive to velocity jumps. A joint inversion with the receiver function allows us to improve theresolution of the seismic discontinuities (e.g., Bodin et al.
2012; Gao and Shen 2015b; Zhou et al. 2020; Akuhara
et al. 2023; Ai et al. 2023). A joint inversion of multimodedispersion with the receiver function for earthquake
data (e.g., Taira and Yoshizawa 2020) could be feasible
to improve depth resolution, although not yet applied to
ANT.

8 Conclusions
The first section introduces theories of SI for ANT.
Although they are based on different assumptions, theresultant CCFs can be related to the Green’s function
under the equipartition of energy. In a realistic situation, although the modal energy is not equally distributed among different mode branches, they work for theindividual mode branches. Heterogeneities of the noise
sources break the equipartition of energy even for a single-mode branch. We evaluated the apparent anomaly
caused by the source heterogeneities. The important
point is that the CCF still satisfies the wave equation evenfor a heterogeneous source distribution, leading to therobustness of the ANT results.
We summarize practical data processing on the calculation of CCFs. For calculating a good quality CCF (satisfying the theoretical assumptions), it is important to choose
appropriate data processing depending on the behavior ofthe data and the noise. Because recent phase velocity measurement methods used waveform information for precise
measurements, amplitude information should be considered by careful data selection and normalization of CCFs.

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Since most conventional ANT studies utilized only thefundamental modes, it has inherited a great uncertainty
regarding the depth structure. Recent developments inmeasurement techniques and dense observations enable us to utilize information on multimode dispersion.
We show a typical inversion of multimode ANT, which
consists of the following four steps. The first step is themultimode dispersion measurement of surface waves
to construct local 1-D structures. The second step is tomeasure the phase/group velocity for each path. Multimode interference was problematic in measuringmultimode dispersion because most studies implicitly
assumed the dominance of the fundamental mode. The
dispersion measurements by waveform fitting in a model
space are feasible for the multimode case. The third stepis the 2-D phase/group velocity inversions. Ray-theoretical inversion was commonly used in most cases, but
the finite-frequency effect has been considered recently.
To evaluate the effects of source heterogeneities on thesensitivity kernel, we show an analytic kernel with source
heterogeneities in a simplified case. Helmholtz tomography is also feasible for considering finite-frequency
effects and the source heterogeneities for sufficiently
dense stations compared to the wavelength because
CCFs for a heterogeneous source distribution still satisfy
the wave equation. The fourth step is a local 1-D inversion on each grid to construct a 3-D S-wave tomographic
model. The multimode inversion improves the depth resolution of the S-wave velocity structure. Joint inversions
with other geophysical observations, such as the H/Vratio and receiver functions, could be feasible to improve
depth resolution further.

Page 35 of 42

MCMC	Markov Chain Monte Carlo
MS	Mean Squared amplitude
multimode ANT	Ambient noise multimode surface wave tomography
NLNM	New Low Noise Model (Peterson 1993)
PREM	
Preliminary Reference Earth Model (Dziewonski andAnderson 1981)
rj-MCMC	Reversible jump Markov Chain Monte Carlo
RMS	Root Mean Squared amplitude
SI	Seismic Interferometry
SPAC	SPatial AutoCorrelation method
USArray	A 15-year program to place a dense network of permanent and portable seismographs across the continental
USA. Seismographs record local, regional, and distant (teleseismic) earthquakes. https://​doi.​org/​10.​7914/​SN/​TA
Acknowledgements
The authors thank many people for maintaining the networks and making thedata readily available. The authors thank N. Nakata, K. Yoshizawa, and T. Tanimoto for constructive comments that significantly improved the quality of themanuscript. The preprocessing and postprocessing of the data analysis werecarried out using ObsPy (Krischer et al. 2015), NumPy (van der Walt et al. 2011),
pyGMT (Tian et al. 2023, Wessel et al. 2019), and SciPy (Virtanen et al. 2020).
Author Contributions
KN contributed to conceptualization; RT and AT contributed to the writing of
the original draft preparation. All authors read, commented on, and approved
the final manuscript.
Funding
K. Nishida was supported by JSPS KAKENHI (Grand No. 21K03715), R. Takagi
was supported by JSPS KAKENHI (Grand No. 23K03522), and A. Takeo was
supported by JSPS KAKENHI (Grand No. 21K14009).
Availability of data and materials
The Python code of the synthetic tests (Figs. 7, 9, 10 and 16) will be available
on the Zenodo page https://​doi.​org/​10.​5281/​zenodo.​10313​732.

Declarations
Conflict of interests
The authors declare that they have no conflict of interest.
Received: 29 September 2023 Accepted: 23 December 2023

Abbreviations
1-D, 2-D, 3-D	One-dimensional, two-dimensional, and therein-dimensional
ANT	Ambient noise surface wave tomography
CCF	
Cross-Correlation Function, which is the time-domain
representation of a cross-spectrum
CMT	Centroid Moment Tensor
FFT	Fast Fourier transform
FJ method	The frequency–Bessel method
FTAN	Frequency–time analysis
Hi-net	
A high-sensitivity seismograph network of nearly 800
stations with an average spacing of 20 km. Seismometers are installed at the bottom of boreholes at a depthof 100–3500 m to reduce noise caused by winds, ocean
waves, and human activity. Observed data is digitizedat each station and continuously transmitted to NIED.
https://​doi.​org/​10.​17598/​NIED.​0003
H/V spectral ratio	The ratio between the Fourier amplitude spectra of thehorizontal and vertical components of microtremors.
IRIS	IRIS (Incorporated Research Institutions for Seismology[1])was a university research consortium dedicated to exploring the interior of the Earth by collecting and distributing
seismographic data. On January 1, 2023, IRIS merged withUNAVCO to form the EarthScope Consortium.
LPO	Lattice-preferred orientationmHz	Milli-Hertz ­(10−3 Hz)

References
Abramowitz M, Stegun IA, Romer RH (1988) Handbook of mathematical
functions with formulas, graphs, and mathematical tables. Am J Phys
56(10):958. https://​doi.​org/​10.​1119/1.​15378
Ai S, Akuhara T, Morishige M, Yoshizawa K, Shinohara M, Nakahigashi K(2023) Layered evolution of the oceanic lithosphere beneath theJapan basin, the sea of Japan. J Geophys Res. https://​doi.​org/​10.​1029/​2022j​b0255​81
Aki K (1955) Relay computer for correlogram analysis. Zisin 8(1):56–58
Aki K (1957) Space and time spectra of stationary stochastic waves, with special reference to microseisms. Bull Earthq Res Inst 35:415–457
Aki K (2003) Earthquake prediction by monitoring the loading process ofplate-driving forces at the Brittle-Ductile transition zone. In: Seventh
workshop on non-linear dynamics and earthquake prediction. Trieste,
Italy, pp 1–166
Aki K, Kaminuma K (1963) Phase velocity of low waves in Japan (part 1):
Love waves from the Aleutian shock of march 9, 1957. PhD thesis, The
University of Tokyo
Aki K, Richards PG (1980) Quantitative seismology, vol 2. San Francisco, W. H.
Freeman, p 557

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Akuhara T (2022) SEISmological transdimensional inversion tools for Flat andIsotropic Layered structures in the Ocean (SEIS\_FILO). https://​doi.​org/​
10.​5281/​zenodo.​63308​40
Akuhara T, Yamashita Y, Ohyanagi S, Sawaki Y, Yamada T, Shinohara M (2023)
Shallow low-velocity layer in the Hyuga-nada accretionary prism andits hydrological implications: insights from a passive seismic array. J
Geophys Res. https://​doi.​org/​10.​1029/​2022j​b0262​98
Anderson DL (1962) Love wave dispersion in heterogeneous anisotropic
media. Geophysics 27(4):445–454. https://​doi.​org/​10.​1190/1.​14390​42
Ardhuin F, Gualtieri L, Stutzmann E (2015) How ocean waves rock the earth:
Two mechanisms explain microseisms with periods 3 to 300 s. Geophys
Res Lett 42(3):765–772. https://​doi.​org/​10.​1002/​2014G​L0627​82
Bahavar M, Spica ZJ, Sánchez-Sesma FJ, Trabant C, Zandieh A, Toro G (2020)
Horizontal-to-vertical spectral ratio (HVSR) IRIS station toolbox. Seismol
Res Lett 91(6):3539–3549. https://​doi.​org/​10.​1785/​02202​00047
Barmin MP (2023) RayTomo: seismic surface wave isotropic/azimuthally anisotropic group/phase speed 2-D tomography. https://​github.​com/​Noise​
CIEI/​RayTo​mo. Accessed 9 Aug 2023
Barmin MP, Ritzwoller MH, Levshin AL (2001) A fast and reliable method forsurface wave tomography. Pure Appl Geophys 158:1351–1375. https://​doi.​org/​10.​1007/​PL000​01225
Barmine M (2023) AFTAN: seismic surface wave dispersion measurementswith automatic frequency time analysis. https://​github.​com/​Noise​CIEI/​AFTAN. Accessed 12 July 2023
Bensen GD, Ritzwoller MH, Barmin MP, Levshin AL, Lin F, Moschetti MP, Shapiro
NM, Yang Y (2007) Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. Geophys J Int
169(3):1239–1260. https://​doi.​org/​10.​1111/j.​1365-​246X.​2007.​03374.x
Bodin T, Sambridge M (2009) Seismic tomography with the reversible jump
algorithm. Geophys J Int 178(3):1411–1436. https://​doi.​org/​10.​1111/j.​
1365-​246X.​2009.​04226.x
Bodin T, Sambridge M, Tkalćić H, Arroucau P, Gallagher K, Rawlinson N (2012)
Transdimensional inversion of receiver functions and surface wave
dispersion. J Geophys Res 117(2):1–24. https://​doi.​org/​10.​1029/​2011J​
B0085​60
Born M, Wolf E, Bhatia AB, Clemmow PC, Gabor D, Stokes AR, Taylor AM,
Wayman PA, Wilcock WL, (1999) Principles of optics: electromagnetic
theory of propagation, interference and diffraction of light, 7th edn.
Cambridge University Press, Cambridge, p 985. https://​doi.​org/​10.​1017/​
CBO97​81139​644181
Boué P, Denolle M, Hirata N, Nakagawa S, Beroza GC (2016) Beyond basin
resonance: characterizing wave propagation using a dense array andthe ambient seismic field. Geophys J Int 206(2):1261–1272. https://​doi.​
org/​10.​1093/​gji/​ggw205
Bowden DC, Tsai VC, Lin FC (2015) Site amplification, attenuation, and scattering from noise correlation amplitudes across a dense array in longbeach. CA Geophys Res Lett 42(5):1360–1367. https://​doi.​org/​10.​1002/​
2014g​l0626​62
Brenguier F, Campillo M, Hadziioannou C, Shapiro NM, Nadeau RM, Larose E(2008) Postseismic relaxation along the San Andreas fault at Parkfieldfrom continuous seismological observations. Science 321(5895):1478–
1481. https://​doi.​org/​10.​1126/​scien​ce.​11609​43
Brenguier F, Shapiro NM, Campillo M, Ferrazzini V, Duputel Z, Coutant O, Nercessian A (2008) Towards forecasting volcanic eruptions using seismic
noise. Nat Geosci 1(2):126–130. https://​doi.​org/​10.​1038/​ngeo1​04
Brocher TM (2005) Empirical relations between elastic wavespeeds and density in the earth’s crust. Bull Seismol Soc Am 95(6):2081–2092. https://​doi.​org/​10.​1785/​01200​50077
Callen HB, Welton TA (1951) Irreversibility and generalized noise. Phys Rev
83(1):34–40. https://​doi.​org/​10.​1103/​PhysR​ev.​83.​34
Campillo M, Paul A (2003) Long-range correlations in the diffuse seismic coda.
Science 299(5606):547–549. https://​doi.​org/​10.​1126/​scien​ce.​10785​51
Chávez-García FJ, Luzón F (2005) On the correlation of seismic microtremors. J
Geophys Res 110(B11):1–12. https://​doi.​org/​10.​1029/​2005J​B0036​71
Chen Y, de Ridder SAL, Rost S, Guo Z, Wu X, Chen Y (2022) Eikonal tomography
with physics-informed neural networks: Rayleigh wave phase velocity
in the northeastern margin of the Tibetan plateau. Geophys Res Lett.
https://​doi.​org/​10.​1029/​2022g​l0990​53
Chmiel M, Roux P, Herrmann P, Rondeleux B, Wathelet M (2018) Data-based
diffraction kernels for surface waves from convolution and correlation

Page 36 of 42

processes through active seismic interferometry. Geophys J Int
214(2):1468–1480. https://​doi.​org/​10.​1093/​gji/​ggy211
Chmiel M, Mordret A, Boué P, Brenguier F, Lecocq T, Courbis R, Hollis D, Campman X, Romijn R, Van der Veen W (2019) Ambient noise multimode
Rayleigh and Love wave tomography to determine the shear velocity
structure above the Groningen gas field. Geophys J Int 218(3):1781–
1795. https://​doi.​org/​10.​1093/​gji/​ggz237
Cho I, Tada T, Shinozaki Y (2006) Centerless circular array method: inferring
phase velocities of Rayleigh waves in broad wavelength ranges using
microtremor records. J Geophys Res. https://​doi.​org/​10.​1029/​2005j​
b0042​35
Claerbout JF (1968) Synthesis of a layered medium from its acoustic transmission response. Geophysics 33(2):264–269. https://​doi.​org/​10.​1190/1.​
14399​27
Cotte N, Laske G (2002) Testing group velocity maps for Eurasia. Geophys J Int
150(3):639–650. https://​doi.​org/​10.​1046/j.​1365-​246X.​2002.​01719.x
Covellone BM, Savage B, Shen Y (2015) Seismic wave speed structure of theOntong Java Plateau. Earth Planet Sci Lett 420:140–150. https://​doi.​org/​
10.​1016/j.​epsl.​2015.​03.​033
Cox H (1973) Spatial correlation in arbitrary noise fields with application toambient sea noise. J Acoust Soc Am 54(5):1289–1301. https://​doi.​org/​
10.​1121/1.​19144​26
Crosbie KJ, Abers GA, Mann ME, Janiszewski HA, Creager KC, Ulberg CW, MoranSC (2019) Shear velocity structure from ambient noise and teleseismic
surface wave tomography in the cascades around mount St Helens. J
Geophys Res 124(8):8358–8375. https://​doi.​org/​10.​1029/​2019j​b0178​36
Cupillard P, Stehly L, Romanowicz B (2011) The one-bit noise correlation:a theory based on the concepts of coherent and incoherent noise.
Geophys J Int 184(3):1397–1414. https://​doi.​org/​10.​1111/j.​1365-​246X.​
2010.​04923.x
Dahlen FA, Tromp J (1998) Theoretical global seismology. Princeton University
Press, Princeton, p 1025
de Ridder S, Dellinger J (2011) Ambient seismic noise Eikonal tomography fornear-surface imaging at Valhall. Lead Edge 30(5):506–512. https://​doi.​
org/​10.​1190/1.​35891​08
de Vos D, Paulssen H, Fichtner A (2013) Finite-frequency sensitivity kernels for
two-station surface wave measurements. Geophys J Int 194(2):1042–
1049. https://​doi.​org/​10.​1093/​gji/​ggt144
Dou S, Lindsey N, Wagner AM, Daley TM, Freifeld B, Robertson M, Peterson
J, Ulrich C, Martin ER, Ajo-Franklin JB (2017) Distributed acoustic
sensing for seismic monitoring of the near surface: a traffic-noiseinterferometry case study. Sci Rep 7(1):11620. https://​doi.​org/​10.​1038/​
s41598-​017-​11986-4
Duvall TL, Jeffferies SM, Harvey JW, Pomerantz MA (1993) Time-distance
helioseismology. Nature 362(6419):430–432. https://​doi.​org/​10.​1038/​
36243​0a0
Dziewonski A, Anderson D (1981) Preliminary reference earth model. Phys
Earth Planet Inter 25:297–356
Dziewonski AM, Hales AL (1972) Numerical analysis of dispersed seismic
waves. In: Bolt BA (ed) Methods in computational physics: advancesin research and applications, vol 11. Elsevier, Amsterdam, pp 39–85.
https://​doi.​org/​10.​1016/​B978-0-​12-​460811-​5.​50007-6
Eckart C (1953) The theory of noise in continuous media. J Acoust Soc Am
25(2):195–199. https://​doi.​org/​10.​1121/1.​19070​18
Ekström G, Tromp J, Larson EWF (1997) Measurements and global models of
surface wave propagation. J Geophys Res 19(17):15705–15710
Ekström G, Ga Abers, Webb SC (2009) Determination of surface-wave phasevelocities across USArray from noise and Aki’s spectral formulation.
Geophys Res Lett 36(18):5–9. https://​doi.​org/​10.​1029/​2009G​L0391​31
Ekström G, Nettles M, Dziewoński AM (2012) The global CMT project 2004–2010: centroid-moment tensors for 13,017 earthquakes. Phys Earth
Planet Inter 200–201:1–9. https://​doi.​org/​10.​1016/j.​pepi.​2012.​04.​002
Emry EL, Shen Y, Nyblade AA, Flinders A, Bao X (2019) Upper mantle earth
structure in Africa from full-wave ambient noise tomography. Geochem
Geophys Geosyst 20(1):120–147. https://​doi.​org/​10.​1029/​2018g​c0078​04
Fichtner A (2014) Source and processing effects on noise correlations. Geophys J Int 197(3):1527–1531. https://​doi.​org/​10.​1093/​gji/​ggu093
Fichtner A (2015) Source-structure trade-offs in ambient noise correlations.
Geophys J Int 202(1):678–694. https://​doi.​org/​10.​1093/​gji/​ggv182

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Fichtner A, Kennett BLN, Trampert J (2013) Separating intrinsic and apparent
anisotropy. Phys Earth Planet Inter 219:11–20. https://​doi.​org/​10.​1016/j.​
pepi.​2013.​03.​006
Fichtner A, Stehly L, Ermert L, Boehm C (2016) Generalised interferometry—I.
Theory for inter-station correlations. Geophys J Int. https://​doi.​org/​10.​1093/​gji/​ggw420
Froment B, Campillo M, Roux P, Gouédard P, Verdel A, Weaver RL (2010)
Estimation of the effect of nonisotropically distributed energy on theapparent arrival time in correlations. Geophysics 75(5):85–93. https://​doi.​org/​10.​1190/1.​34831​02
Fu L, Pan L, Li Z, Dong S, Ma Q, Chen X (2022) Improved high-resolution 3D vsmodel of long beach, CA: inversion of multimodal dispersion curves
from ambient noise of a dense array. Geophys Res Lett. https://​doi.​org/​
10.​1029/​2021g​l0976​19
Fukao Y, Nishida K, Suda N, Nawa K, Kobayashi N (2002) A theory of the earth’s
background free oscillations. J Geophys Res 107(B9):2206
Fukao Y, Nishida K, Kobayashi N (2010) Seafloor topography, ocean infragravity waves, and background Love and Rayleigh waves. J Geophys Res
115(B4):04302. https://​doi.​org/​10.​1029/​2009J​B0066​78
Fukushima S, Shinohara M, Nishida K, Takeo A, Yamada T, Yomogida K (2022)
Detailed s-wave velocity structure of sediment and crust off Sanriku,
Japan by a new analysis method for distributed acoustic sensing datausing a seafloor cable and seismic interferometry. Earth Planets Space
74(1):1–11. https://​doi.​org/​10.​1186/​s40623-​022-​01652-z
Gabriels P, Snieder R, Nolet G (1987) In situ measurements of shear-wave
velocity in sediments with higher-mode Rayleigh waves. Geophys Prospect 35(2):187–196. https://​doi.​org/​10.​1111/j.​1365-​2478.​1987.​tb008​12.x
Galetti E, Curtis A, Baptie B, Jenkins D, Nicolson H (2016) Transdimensional
Love-wave tomography of the British isles and shear-velocity structureof the east Irish sea basin from ambient-noise interferometry. Geophys
J Int. https://​doi.​org/​10.​1093/​gji/​ggw286
Gao H, Shen Y (2014) Upper mantle structure of the cascades from full-wave
ambient noise tomography: evidence for 3D mantle upwelling in the
back-arc. Earth Planet Sci Lett 390:222–233. https://​doi.​org/​10.​1016/j.​
epsl.​2014.​01.​012
Gao H, Shen Y (2015) A preliminary Full-Wave Ambient-Noise tomographymodel spanning from the Juan de Fuca and Gorda spreading centers to
the Cascadia volcanic arc. Seismol Res Lett 86(5):1253–1260. https://​doi.​
org/​10.​1785/​02201​50103
Gao H, Shen Y (2015) Validation of recent shear wave velocity models in the
United States with full-wave simulation. J Geophys Res 120(1):344–358.
https://​doi.​org/​10.​1002/​2014j​b0113​69
Gizon L, Birch AC (2002) Time–distance helioseismology: the forward problem
for random distributed sources. Astrophys J 571(2):966–986. https://​doi.​
org/​10.​1086/​340015
Godin OA (2007) Emergence of the acoustic Green’s function from thermal
noise. J Acoust Soc Am 121(2):96–102. https://​doi.​org/​10.​1121/1.​24307​
64
Gouédard P, Stehly L, Brenguier F, Campillo M, Colin de Verdière Y, Larose E,
Margerin L, Roux P, Sánchez-Sesma FJ, Shapiro NM, Weaver RL (2008)
Cross-correlation of random fields: mathematical approach and applications. Geophys Prospect 56(3):375–393. https://​doi.​org/​10.​1111/j.​
1365-​2478.​2007.​00684.x
Gouédard P, Seher T, McGuire JJ, Collins JA, van der Hilst RD (2014) Correctionof ocean-bottom seismometer instrumental clock errors using ambient
seismic noise. Bull Seism Soc Am 104(3):1276–1288
Green PJ (1995) Reversible jump Markov chain Monte Carlo computation andBayesian model determination. Biometrika 82(4):711–732. https://​doi.​
org/​10.​1093/​biomet/​82.4.​711
Hamada K, Yoshizawa K (2015) Interstation phase speed and amplitude measurements of surface waves with nonlinear waveform fitting: application
to USArray. Geophys J Int 202(3):1463–1482. https://​doi.​org/​10.​1093/​
gji/​ggv213
Hanasoge SM (2013) Measurements and kernels for source-structure inversions in noise tomography. Geophys J Int 196(2):971–985. https://​doi.​
org/​10.​1093/​gji/​ggt411
Hanasoge S, Gizon L, Sreenivasan KR (2016) Seismic sounding of convection in
the sun. Annu Rev Fluid Mech 48(1):191–217. https://​doi.​org/​10.​1146/​
annur​ev-​fluid-​122414-​034534

Page 37 of 42

Haned A, Stutzmann E, Schimmel M, Kiselev S, Davaille A, Yelles-Chaouche
A (2016) Global tomography using seismic hum. Geophys J Int
204(2):1222–1236. https://​doi.​org/​10.​1093/​gji/​ggv516
Haney MM, Mikesell TD, van Wijk K, Nakahara H (2012) Extension of the spatial
autocorrelation (SPAC) method to mixed-component correlations ofsurface waves. Geophys J Int 191(1):189–206. https://​doi.​org/​10.​1111/j.​
1365-​246X.​2012.​05597.x
Hannemann K, Krüger F, Dahm T (2014) Measuring of clock drift rates andstatic time offsets of ocean bottom stations by means of ambient noise.
Geophys J Int 196(2):1034–1042
Harmon N, Rychert CA (2015) Seismic imaging of deep crustal melt sills
beneath Costa Rica suggests a method for the formation of theArchean continental crust. Earth Planet Sci Lett 430:140–148. https://​doi.​org/​10.​1016/j.​epsl.​2015.​07.​062
Harmon N, Forsyth D, Webb S (2007) Using ambient seismic noise todetermine short-period phase velocities and shallow shear velocities
in young oceanic lithosphere. Bull Seismol Soc Am 97(6):2009–2023.
https://​doi.​org/​10.​1785/​01200​70050
Harmon N, Rychert C, Gerstoft P (2010) Distribution of noise sources for seismic interferometry. Geophys J Int 183(3):1470–1484. https://​doi.​org/​10.​1111/j.​1365-​246X.​2010.​04802.x
Hasselmann K (1963) A statistical analysis of the generation of microseisms.
Rev Geophys 1(2):177. https://​doi.​org/​10.​1029/​RG001​i002p​00177
Herrmann RB (2013) Computer programs in seismology: an evolving tool for
instruction and research. Seismol Res Lett 84(6):1081–1088. https://​doi.​
org/​10.​1785/​02201​10096
Hirose T, Ueda H (2023) Continuous monitoring of instrumental clock errors at
50 volcanoes in Japan based on seismic interferometry. Earth Planets
Space 75(1):1–12. https://​doi.​org/​10.​1186/​s40623-​023-​01798-4
Hirose T, Nakahara H, Nishimura T, Campillo M (2020) Locating spatial changesof seismic scattering property by sparse modeling of seismic ambient
noise cross-correlation functions: application to the 2008 Iwate-Miyagi
Nairiku (M w 6.9), Japan, earthquake. J Geophys Res. https://​doi.​org/​10.​1029/​2019j​b0193​07
Hu S, Luo S, Yao H (2020) The frequency-Bessel spectrograms of multicomponent cross-correlation functions from seismic ambient noise. J Geophys
Res. https://​doi.​org/​10.​1029/​2020j​b0196​30
Huang H, Yao H, van der Hilst RD (2010) Radial anisotropy in the crust of SE
Tibet and SW China from ambient noise interferometry. Geophys Res
Lett 37(21):1–5. https://​doi.​org/​10.​1029/​2010G​L0449​81
Ikeda T, Matsuoka T, Tsuji T, Hayashi K (2012) Multimode inversion with amplitude response of surface waves in the spatial autocorrelation method.
Geophys J Int 190(1):541–552. https://​doi.​org/​10.​1111/j.​1365-​246X.​
2012.​05496.x
IRIS DMC Data Services Products: EMC, A repository of Earth model. https://​doi.​org/​10.​17611/​DP/​EMC.1
Ishimaru A (1997) Wave propagation and scattering in random media. IEEE
Press, New York
Janiszewski HA, Gaherty JB, Abers GA, Gao H, Eilon ZC (2019) Amphibious
surface-wave phase-velocity measurements of the Cascadia subduction zone. Geophys J Int 217(3):1929–1948. https://​doi.​org/​10.​1093/​gji/​ggz051
Jaxybulatov K, Shapiro NM, Koulakov I, Mordret A, Landes M, Sens-Schonfelder
C (2014) A large magmatic sill complex beneath the Toba caldera. Science 346(6209):617–619. https://​doi.​org/​10.​1126/​scien​ce.​12585​82
Jiang C, Denolle MA (2022) Pronounced seismic anisotropy in Kanto sedimentary basin: a case study of using dense arrays, ambient noise seismology, and multi-modal surface-wave imaging. J Geophys Res. https://​doi.​org/​10.​1029/​2022j​b0246​13
Jiang C, Schmandt B, Farrell J, Lin F-C, Ward KM (2018) Seismically anisotropic
magma reservoirs underlying silicic calderas. Geology 46(8):727–730.
https://​doi.​org/​10.​1130/​G45104.1
Jiang C, Schmandt B, Abers GA, Kiser E, Miller MS (2023) Segmentation andradial anisotropy of the deep crustal magmatic system beneath thecascades arc. Geochem Geophys Geosyst. https://​doi.​org/​10.​1029/​
2022g​c0107​38
Kawano Y, Isse T, Takeo A, Kawakatsu H, Morishige M, Shiobara H, Takeuchi
N, Sugioka H, Kim Y, Utada H, Lee S-M (2023) Seismic structure of
the lithosphere-asthenosphere system beneath the oldest seafloor

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

revealed by Rayleigh-wave dispersion analysis. J Geophys Res. https://​doi.​org/​10.​1029/​2023j​b0265​29
Kobayashi N, Nishida K (1998) Continuous excitation of planetary free oscillations by atmospheric disturbances. Nature 395(September):357–360.
https://​doi.​org/​10.​1038/​26427
Krischer L, Megies T, Barsch R, Beyreuther M, Lecocq T, Caudron C, Wassermann
J (2015) ObsPy: a bridge for seismology into the scientific python
ecosystem. Comput Sci Discov 8(1):014003. https://​doi.​org/​10.​1088/​
1749-​4699/8/​1/​014003
Kubo R (1957) Statistical-Mechanical theory of irreversible processes. I. general
theory and simple applications to magnetic and conduction problems.
J Phys Soc Jpn 12(6):570–586. https://​doi.​org/​10.​1143/​JPSJ.​12.​570
Landisman M, Dziewonski A, Satô Y (1969) Recent improvements in the analysis of surface wave observations. Geophys J Int 17(4):369–403. https://​doi.​org/​10.​1111/j.​1365-​246X.​1969.​tb002​46.x
Laske G, Widmer-Schnidrig R (2015) 1.04—theory and observations: normal
mode and surface wave observations. In: Schubert G (ed) Treatise on
geophysics, 2nd edn. Elsevier, Oxford, pp 117–167. https://​doi.​org/​10.​1016/​B978-0-​444-​53802-4.​00003-8
Lehujeur M, Chevrot S (2020) Eikonal tomography using coherent surface
waves extracted from ambient noise by iterative matched filtering—
application to the large-N Maupasacq array. J Geophys Res. https://​doi.​
org/​10.​1029/​2020j​b0193​63
Levshin AL, Barmin MP, Ritzwoller MH (2018) Tutorial review of seismic surface
waves’ phenomenology. J Seismol 22(2):519–537. https://​doi.​org/​10.​
1007/​s10950-​017-​9716-7
Levshin, A SensKernel: 1-D seismic surface wave isotropicgroup/phase speed
dispersion/eigenfunction/sensitivity kernel. https://​github.​com/​Noise​
CIEI/​SensK​ernel
Levshin A, Ratnikova L, Berger J (1992) Peculiarities of surface-wave propagation across central Eurasia. Bull Seismol Soc Am 82(6):2464–2493.
https://​doi.​org/​10.​1785/​BSSA0​82006​2464
Li J, Chen Y, Schuster GT (2020) Separation of multi-mode surface waves bysupervised machine learning methods. Geophys Prospect 68(4):1270–
1280. https://​doi.​org/​10.​1111/​1365-​2478.​12927
Liang C, Langston C (2008) Ambient seismic noise tomography and structureof eastern North America. J Geophys Res 113(B3):1–18. https://​doi.​org/​
10.​1029/​2007J​B0053​50
Lin F-C, Ritzwoller MH (2011) Apparent anisotropy in inhomogeneous
isotropic media. Geophys J Int 186(3):1205–1219. https://​doi.​org/​10.​1111/j.​1365-​246X.​2011.​05100.x
Lin FC, Ritzwoller MH (2011) Helmholtz surface wave tomography for isotropic
and azimuthally anisotropic structure. Geophys J Int 186(3):1104–1120.
https://​doi.​org/​10.​1111/j.​1365-​246X.​2011.​05070.x
Lin F-C, Schmandt B (2014) Upper crustal azimuthal anisotropy across thecontiguous U.S. determined by Rayleigh wave ellipticity. Geophys Res
Lett 41(23):8301–8307. https://​doi.​org/​10.​1002/​2014g​l0623​62
Lin F-C, Moschetti MP, Ritzwoller MH (2008) Surface wave tomography of thewestern United States from ambient seismic noise: Rayleigh and Lovewave phase velocity maps. Geophys J Int 173(1):281–298. https://​doi.​
org/​10.​1111/j.​1365-​246X.​2008.​03720.x
Lin F-C, Ritzwoller MH, Snieder R (2009) Eikonal tomography: surface wave
tomography by phase front tracking across a regional broad-bandseismic array. Geophys J Int 177(3):1091–1110. https://​doi.​org/​10.​1111/j.​
1365-​246X.​2009.​04105.x
Lin F-C, Ritzwoller MH, Shen W (2011) On the reliability of attenuation measurements from ambient noise cross-correlations. Geophys Res Lett.
https://​doi.​org/​10.​1029/​2011g​l0473​66
Lin F-C, Schmandt B, Tsai VC (2012) Joint inversion of Rayleigh wave phasevelocity and ellipticity using USArray: constraining velocity and density
structure in the upper crust. Geophys Res Lett 39:1–7. https://​doi.​org/​
10.​1029/​2012G​L0521​96
Lin F-C, Tsai VC, Ritzwoller MH (2012) The local amplification of surface waves:a new observable to constrain elastic velocities, density, and anelastic
attenuation. J Geophys Res 117:B06302. https://​doi.​org/​10.​1029/​2012J​
B0092​08
Lin F-C, Li D, Clayton RW, Hollis D (2013) High-resolution 3D shallow crustal
structure in long beach, California: application of ambient noise tomography on a dense seismic array. Geophysics 78(4):45–56. https://​doi.​org/​
10.​1190/​geo20​12-​0453.1

Page 38 of 42

Lin P-YP, Gaherty JB, Jin G, Collins JA, Lizarralde D, Evans RL, Hirth G (2016)
High-resolution seismic constraints on flow dynamics in the oceanic
asthenosphere. Nature 535(7613):538–541. https://​doi.​org/​10.​1038/​
natur​e18012
Liu X (2020) Finite-frequency sensitivity kernels for seismic noise interferometry based on differential time measurements. J Geophys Res. https://​doi.​org/​10.​1029/​2019j​b0189​32
Liu X, Ben-Zion Y (2013) Theoretical and numerical results on effects of attenuation on correlation functions of ambient seismic noise. Geophys J Int
194(3):1966–1983. https://​doi.​org/​10.​1093/​gji/​ggt215
Liu X, Beroza GC, Yang L, Ellsworth WL (2021) Ambient noise Love wave
attenuation tomography for the LASSIE array across the Los Angeles
basin. Sci Adv. https://​doi.​org/​10.​1126/​sciadv.​abe10​30
Liu Q, Chen X, Gao L, Yu Z, Chen J (2023) Direct image dissimilarity inversion ofambient noise multimodal dispersion spectrograms. Bull Seismol Soc
Am. https://​doi.​org/​10.​1785/​01202​30012
Lobkis OI, Weaver RL (2001) On the emergence of the green’s function in thecorrelations of a diffuse field. J Acoust Soc Am 110(6):3011. https://​doi.​
org/​10.​1121/1.​14175​28
Longuet-Higgins MS (1950) A theory of the origin of microseisms. Philos TransR Soc Lond 243(857):1–35. https://​doi.​org/​10.​1098/​rsta.​1950.​0012
Lu, Y, Stehly, L, Brossier, R, Paul, A, AlpArray Working Group (2020) Imaging
alpine crust using ambient noise wave-equation tomography. Geophys
J Int 222(1):69–85. https://​doi.​org/​10.​1093/​gji/​ggaa1​45
Lynner C, Beck SL, Zandt G, Porritt RW, Lin FC, Eilon ZC (2018) Midcrustal
deformation in the central Andes constrained by radial anisotropy. J
Geophys Res 123(6):4798–4813. https://​doi.​org/​10.​1029/​2017J​B0149​36
Magrini F, Boschi L (2021) Surface-wave attenuation from seismic ambient
noise: numerical validation and application. J Geophys Res. https://​doi.​
org/​10.​1029/​2020j​b0198​65
Maraschini M, Foti S (2010) A Monte Carlo multimodal inversion of surface
waves. Geophys J Int 182(3):1557–1566. https://​doi.​org/​10.​1111/j.​1365-​
246X.​2010.​04703.x
Menke W, Jin G (2015) Waveform fitting of cross spectra to determine phasevelocity using Aki’s formula. Bull Seismol Soc Am 105(3):1619–1627.
https://​doi.​org/​10.​1785/​01201​40245
Miller D, Bennington N, Haney M, Bedrosian P, Key K, Thurber C, Hart L, Ohlendorf S (2020) Linking magma storage and ascent to eruption volume
and composition at an arc caldera. Geophys Res Lett. https://​doi.​org/​10.​1029/​2020g​l0881​22
Montagner J-P (2015) Deep earth structure—upper mantle structure: global
isotropic and anisotropic elastic tomography. In: Schubert G (ed) Treatise on geophysics, vol 1, 2nd edn. Elsevier, Amsterdam, pp 613–639.
https://​doi.​org/​10.​1016/​B978-0-​444-​53802-4.​00016-6
Mordret A, Shapiro NM, Singh SS, Roux P, Barkved OI (2013) Helmholtz
tomography of ambient noise surface wave data to estimate Scholte
wave phase velocity at Valhall life of the field. Geophysics 78(2):99–109.
https://​doi.​org/​10.​1190/​geo20​12-​0303.1
Moschetti MP, Ritzwoller MH, Shapiro NM (2007) Surface wave tomography of
the western United States from ambient seismic noise: Rayleigh wave
group velocity maps. Geochem Geophys Geosyst. https://​doi.​org/​10.​1029/​2007g​c0016​55
Moschetti MP, Ritzwoller MH, Lin F-C, Yang Y (2010) Crustal shear wave velocity
structure of the western United States inferred from ambient seismic
noise and earthquake data. J Geophys Res. https://​doi.​org/​10.​1029/​
2010j​b0074​48
Moschetti MP, Ritzwoller MH, Lin F, Yang Y (2010) Seismic evidence for widespread western-US deep-crustal deformation caused by extension.
Nature 464(7290):885–889. https://​doi.​org/​10.​1038/​natur​e08951
Nagaoka Y, Nishida K, Aoki Y, Takeo M, Ohminato T (2012) Seismic imaging
of magma chamber beneath an active volcano. Earth Planet Sci Lett
333–334:1–8. https://​doi.​org/​10.​1016/j.​epsl.​2012.​03.​034
Nakahara H (2006) A systematic study of theoretical relations between spatial
correlation and Green’s function in one-, two- and three-dimensionalrandom scalar wavefields. Geophys J Int 165(3):1097–1105. https://​doi.​
org/​10.​1111/j.​1365-​246X.​2006.​03170.x
Nakahara H, Emoto K, Nishimura T (2021) Extending the formulation of thespatial autocorrelation (SPAC) method to strain, rotation and tilt. Geophys J Int 227(1):287–302. https://​doi.​org/​10.​1093/​gji/​ggab2​17

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Nakamura Y (1989) A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface. Q Rep RTRI 30(1)Nakanishi I, Anderson DL (1982) Worldwide distribution of group velocity ofmantle Rayleigh waves as determined by spherical harmonic inversion.
Bull Seismol Soc Am 72(4):1185–1194. https://​doi.​org/​10.​1785/​BSSA0​
72004​1185
Nakata N (2020) Pseudo arrivals generated by frequency normalization forseismic interferometry with scattered waves—stationary-phase analysis. SEG Tech Program Expanded Abstr 2020:2085–2089. https://​doi.​org/​
10.​1190/​segam​2020-​34279​03.1
Nakata N, Snieder R, Kuroda S, Ito S, Aizawa T, Kunimi T (2013) Monitoringa building using deconvolution interferometry. I: earthquake-data
analysis. Bull Seismol Soc Am 103(3):1662–1678. https://​doi.​org/​10.​
1785/​01201​20291
Nakata N, Gualtieri L, Fichtner A (2019) Seismic ambient noise. Cambridge
University Press, Cambridge
Nayak A, Thurber CH (2020) Using multicomponent ambient seismic noisecross-correlations to identify higher mode Rayleigh waves and improve
dispersion measurements. Geophys J Int 222(3):1590–1605. https://​doi.​
org/​10.​1093/​gji/​ggaa2​70
Nimiya H, Ikeda T, Tsuji T (2020) Three-dimensional S wave velocity structure ofcentral Japan estimated by surface-wave tomography using ambient
noise. J Geophys Res 125(4):2019–019043. https://​doi.​org/​10.​1029/​
2019j​b0190​43
Nimiya H, Ikeda T, Tsuji T (2023) Multimodal Rayleigh and Love wave joint
inversion for s-wave velocity structures in Kanto basin, Japan. J Geophys
Res. https://​doi.​org/​10.​1029/​2022j​b0250​17
Nishida K (2011) Two-dimensional sensitivity kernels for cross-correlation
functions of background surface waves. C R Geosci 343(8–9):584–590.
https://​doi.​org/​10.​1016/j.​crte.​2011.​02.​004
Nishida K (2013) Earth’s background free oscillations. Annu Rev Earth Planet Sci
41(1):719–740. https://​doi.​org/​10.​1146/​annur​ev-​earth-​050212-​124020
Nishida K (2014) Source spectra of seismic hum. Geophys J Int 199(1):416–429.
https://​doi.​org/​10.​1093/​gji/​ggu272
Nishida K (2017) Ambient seismic wave field. Proc Jpn Acad Ser B Phys Biol Sci
93(7):423–448. https://​doi.​org/​10.​2183/​pjab.​93.​026
Nishida K, Kobayashi N (1999) Statistical features of earth’s continuous freeoscillations. J Geophys Res 104(B12):28741–28
Nishida K, Takagi R (2022) A global centroid single force catalog of p-wave
microseisms. J Geophys Res. https://​doi.​org/​10.​1029/​2021j​b0234​84
Nishida K, Kobayashi N, Fukao Y (2002) Origin of earth’s ground noise from 2 to
20 mHz. Geophys Res Lett 29(10):1413
Nishida K, Kawakatsu H, Obara K (2008) Three-dimensional crustal S wave
velocity structure in japan using microseismic data recorded by hi-net
tiltmeters. J Geophys Res 113(B10):10302. https://​doi.​org/​10.​1029/​2007J​
B0053​95
Nishida K, Kawakatsu H, Fukao Y, Obara K (2008) Background Love andRayleigh waves simultaneously generated at the pacific ocean floors.
Geophys Res Lett 35(16):16307. https://​doi.​org/​10.​1029/​2008g​l0347​53
Nishida K, Montagner J-P, Kawakatsu H (2009) Global surface wave tomography using seismic hum. Science 326(5949):112. https://​doi.​org/​10.​1126/​
scien​ce.​11763​89
Nolet G (2008) A breviary of seismic tomography: imaging the interior of theearth and sun. Cambridge University Press, Cambridge. https://​doi.​org/​
10.​1017/​CBO97​80511​984709
Obermann A, Hillers G (2019) Chapter two—seismic time-lapse interferometry
across scales. In: Schmelzbach C (ed) Advances in geophysics, vol 60. Elsevier, Amsterdam, pp 65–143. https://​doi.​org/​10.​1016/​bs.​agph.​2019.​06.​001
Okada H, Sakajiri N (1983) Estimates of an s-wave velocity distribution using
long-period microtremors. Geophys Bull Hokkaido Univ 42:119–143
Oppenheim AV, Schafer RW (2014) Discrete-time signal processing. Pearson, Essex
Pedersen HA, Mattern F, Poli P, Stehly L (2023) Imaging with seismic noise:
improving extraction of body wave phases from the deep earth
through selective stacking based on H/V ratios. Geophys J Int. https://​doi.​org/​10.​1093/​gji/​ggac3​88
Peterson JR (1993) Observations and modeling of seismic background noise.
U.S. Geol Surv Open File Rep 93–322:1–94. https://​doi.​org/​10.​3133/​
ofr93​322
Prieto GA, Lawrence JF, Beroza GC (2009) Anelastic earth structure from thecoherency of the ambient seismic field. J Geophys Res 114(B7):1–15.
https://​doi.​org/​10.​1029/​2008J​B0060​67

Page 39 of 42

Prieto GA, Denolle M, Lawrence JF, Beroza GC (2011) On amplitude information carried by the ambient seismic field. C R Geosci 343(8–9):600–614.
https://​doi.​org/​10.​1016/j.​crte.​2011.​03.​006
Qin T, Lu L, Ding Z, Feng X, Zhang Y (2022) High-resolution 3D shallow S wave
velocity structure of Tongzhou, subcenter of Beijing, inferred from multimode Rayleigh waves by beamforming seismic noise at a dense array.
J Geophys Res. https://​doi.​org/​10.​1029/​2021j​b0236​89
Rawlinson N (2023) Surface wave tomography code. http://​rses.​anu.​edu.​au/​~nick/​surft​omo.​html. Accessed 9 Aug 2023
Rawlinson N, Sambridge M (2004) Wave front evolution in strongly heterogeneous layered media using the fast marching method. Geophys J Int
156(3):631–647. https://​doi.​org/​10.​1111/j.​1365-​246X.​2004.​02153.x
Rawlinson N, Pozgay S, Fishwick S (2010) Seismic tomography: a window into
deep earth. Phys Earth Planet Inter 178(3–4):101–135. https://​doi.​org/​
10.​1016/j.​pepi.​2009.​10.​002
Ritzwoller MH (2002) Global surface wave diffraction tomography. J Geophys
Res 107(B12):1–13. https://​doi.​org/​10.​1029/​2002J​B0017​77
Ritzwoller M, Feng L (2019) Overview of pre- and post-processing of
ambient-noise correlations. In: Nakata N, Gualtieri L, Fichtner A (eds)
Seismic ambient noise. Cambridge University Press, Cambridge, pp
144–187. https://​doi.​org/​10.​1017/​97811​08264​808.​007
Ritzwoller MH, Lin F-CC, Shen W (2011) Ambient noise tomography with alarge seismic array. C R Geosci 343(8–9):558–570. https://​doi.​org/​10.​1016/j.​crte.​2011.​03.​007
Romanowicz B (2003) Global mantle tomography: progress status in thepast 10 years. Annu Rev Earth Planet Sci 31(1):303–328. https://​doi.​
org/​10.​1146/​annur​ev.​earth.​31.​091602.​113555
Romanowicz B (2020) Surface waves. In: Gupta HK (ed) Encyclopedia ofsolid earth geophysics. Springer, Cham, pp 1–13. https://​doi.​org/​10.​
1007/​978-3-​030-​10475-7_​143-1
Romanowicz B (2021) Seismic tomography of the earth’s mantle. In: Alderton D, Elias SA (eds) Encyclopedia of geology, 2nd edn. Academic
Press, Oxford, pp 587–609. https://​doi.​org/​10.​1016/​B978-0-​08-​
102908-​4.​00169-7
Rost S, Thomas C (2002) Array seismology: methods and applications. Rev
Geophys. https://​doi.​org/​10.​1029/​2000R​G0001​00
Russell JB, Gaherty JB (2021) Lithosphere structure and seismic anisotropy
offshore eastern North America: implications for continental breakup
and ultra-slow spreading dynamics. J Geophys Res. https://​doi.​org/​
10.​1029/​2021j​b0229​55
Sager K, Ermert L, Boehm C, Fichtner A (2017) Towards full waveform ambient noise inversion. Geophys J Int 212(1):566–590. https://​doi.​org/​10.​1093/​gji/​ggx429
Sager K, Boehm C, Ermert L, Krischer L, Fichtner A (2020) Global-scale fullwaveform ambient noise inversion. J Geophys Res. https://​doi.​org/​10.​1029/​2019j​b0186​44
Saito M (1988) DISPER80; a subroutine package for calculation of seismic
normal-mode solutions. Seismological algorithms; computational
methods and computer programs. Academic Press, San Diego, pp
293–319
Sánchez-Sesma FJ, Rodríguez M, Iturrarán-Viveros U, Luzón F, Campillo M,
Margerin L, García-Jerez A, Suarez M, Santoyo MA, Rodríguez-Castellanos A (2011) A theory for microtremor H/V spectral ratio: application for a layered medium. Geophys J Int 186(1):221–225. https://​doi.​
org/​10.​1111/j.​1365-​246X.​2011.​05064.x
Sato H, Fehler MC, Maeda T (2012) Seismic wave propagation and scatteringin the heterogeneous earth. Springer, Heidelberg, p 496
Sato H (2019) Isotropic scattering coefficient of the solid earth. Geophys J
Int 218(3):2079–2088. https://​doi.​org/​10.​1093/​gji/​ggz266
Savage MK, Lin F-C, Townend J (2013) Ambient noise cross-correlation
observations of fundamental and higher-mode Rayleigh wave
propagation governed by basement resonance. Geophys Res Lett
40(14):3556–3561. https://​doi.​org/​10.​1002/​grl.​50678
Saygin E, Kennett BLN (2012) Crustal structure of Australia from ambient
seismic noise tomography. J Geophys Res 117(1):1–15. https://​doi.​
org/​10.​1029/​2011J​B0084​03
Saygin E, Kennett BLN (2010) Ambient seismic noise tomography of Australian continent. Tectonophysics 481(1):116–125. https://​doi.​org/​10.​1016/j.​tecto.​2008.​11.​013
Schuster GT (2009) Seismic interferometry. Cambridge University Press,
Cambridge. https://​doi.​org/​10.​1017/​cbo97​80511​581557

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Sens-Schönfelder C (2008) Synchronizing seismic networks with ambient
noise. Geophys J Int 174(3):966–970. https://​doi.​org/​10.​1111/j.​1365-​
246X.​2008.​03842.x
Sens-Schönfelder C, Wegler U (2006) Passive image interferometry andseasonal variations of seismic velocities at Merapi volcano, Indonesia.
Geophys Res Lett 33(21):21302. https://​doi.​org/​10.​1029/​2006G​L0277​
97
Sens-Schönfelder C, Wegler U (2006) Radiative transfer theory for estimationof the seismic moment. Geophys J Int 167(3):1363–1372. https://​doi.​
org/​10.​1111/j.​1365-​246X.​2006.​03139.x
Sens-Schönfelder C, Wegler U (2011) Passive image interferometry formonitoring crustal changes with ambient seismic noise. C R Geosci343(8–9):639–651. https://​doi.​org/​10.​1016/j.​crte.​2011.​02.​005
Sethian JA (1996) A fast marching level set method for monotonicallyadvancing fronts. Proc Natl Acad Sci U S A 93(4):1591–1595. https://​doi.​org/​10.​1073/​pnas.​93.4.​1591
Shapiro NM, Campillo M (2004) Emergence of broadband Rayleigh waves fromcorrelations of the ambient seismic noise. Geophys Res Lett 31(7):1615–
1619. https://​doi.​org/​10.​1029/​2004G​L0194​91
Shapiro NM, Singh SK (1999) A systematic error in estimating surface-wave
group-velocity dispersion curves and a procedure for its correction.
Bull Seismol Soc Am 89(4):1138–1142. https://​doi.​org/​10.​1785/​BSSA0​
89004​1138
Shapiro NMNM, Ritzwoller MHH (2002) Monte-Carlo inversion for a global
shear-velocity model of the crust and upper mantle. Geophys J Int
151(1):88–105. https://​doi.​org/​10.​1046/j.​1365-​246X.​2002.​01742.x
Shapiro NM, Campillo M, Stehly L, Ritzwoller MH (2005) High-resolution
surface-wave tomography from ambient seismic noise. Science
307(5715):1615–1618. https://​doi.​org/​10.​1126/​scien​ce.​11083​39
Snieder R (2004) Extracting the Green’s function from the correlation of coda
waves: a derivation based on stationary phase. Phys Rev E 69(4):046610.
https://​doi.​org/​10.​1103/​PhysR​evE.​69.​046610
Snieder R (2006) Retrieving the Green’s function of the diffusion equation fromthe response to a random forcing. Phys Rev E 74(4):046620. https://​doi.​
org/​10.​1103/​PhysR​evE.​74.​046620
Snieder R, Larose E (2013) Extracting earth’s elastic wave response from noise
measurements. Annu Rev Earth Planet Sci 41(1):183–206. https://​doi.​
org/​10.​1146/​annur​ev-​earth-​050212-​123936
Snieder R, Wapenaar K (2010) Imaging with ambient noise. Phys Today
63(9):44–49. https://​doi.​org/​10.​1063/1.​34905​00
Snieder R, Sheiman J, Calvert R (2006) Equivalence of the virtual-source
method and wave-field deconvolution in seismic interferometry. Phys
Rev E 73(6 Pt 2):066620. https://​doi.​org/​10.​1103/​PhysR​evE.​73.​066620
Socco LV, Foti S, Boiero D (2010) Surface-wave analysis for building near-surface velocity models—established approaches and new perspectives.
Geophysics 75(5):75–8375102. https://​doi.​org/​10.​1190/1.​34794​91
Spetzler J, Trampert J, Snieder R (2002) The effect of scattering in surface wave
tomography. Geophys J Int 149(3):755–767
Spica Z, Perton M, Legrand D (2017) Anatomy of the Colima volcano magmatic system, Mexico. Earth Planet Sci Lett 459:1–13. https://​doi.​org/​10.​1016/j.​epsl.​2016.​11.​010
Spica Z, Perton M, Nakata N, Liu X, Beroza GC (2018) Shallow VS imaging of
the Groningen area from joint inversion of multimode surface waves
and H/V spectral ratios. Seismol Res Lett. https://​doi.​org/​10.​1785/​02201​
80060
Taira T, Yoshizawa K (2020) Upper-mantle discontinuities beneath Australia
from transdimensional Bayesian inversions using multimode surface
waves and receiver functions. Geophys J Int 223(3):2085–2100. https://​doi.​org/​10.​1093/​gji/​ggaa4​42
Takagi R, Nishida K (2022) Multi-mode dispersion measurement of surface
waves extracted by multi-component ambient noise cross-correlation
functions. Geophys J Int. https://​doi.​org/​10.​1093/​gji/​ggac2​25
Takagi R, Nishida K, Aoki Y, Maeda T, Masuda K, Takeo M, Obara K, Shiomi K,
Sato M, Saito K (2015) A single bit matters: coherent noise of seismic
data loggers. Seismol Res Lett 86(3):901–907. https://​doi.​org/​10.​1785/​
02201​50030
Takagi R, Nakahara H, Kono T, Okada T (2014) Separating body and Rayleighwaves with cross terms of the cross-correlation tensor of ambient
noise. J Geophys Res 119(3):2005–2018. https://​doi.​org/​10.​1002/​2013J​
B0108​24

Page 40 of 42

Takagi R, Nishida K, Maeda T, Obara K (2018) Ambient seismic noise wavefieldin japan characterized by polarization analysis of hi-net records. Geophys J Int 215(3):1682–1699. https://​doi.​org/​10.​1093/​gji/​ggy334
Takagi R, Uchida N, Nakayama T, Azuma R, Ishigami A, Okada T, Nakamura
T, Shiomi K (2019) Estimation of the orientations of the s-net cabled
Ocean-Bottom sensors. Seismol Res Lett 90(6):2175–2187. https://​doi.​
org/​10.​1785/​02201​90093
Takagi R, Toyokuni G, Chikasada N (2021) Ambient noise correlation analysis ofs-net records: extracting surface wave signals below instrument noiselevels. Geophys J Int 224(3):1640–1657. https://​doi.​org/​10.​1093/​gji/​ggaa5​48
Takeo A, Nishida K, Isse T, Kawakatsu H, Shiobara H, Sugioka H, Kanazawa T(2013) Radially anisotropic structure beneath the Shikoku basin frombroadband surface wave analysis of ocean bottom seismometerrecords. J Geophys Res 118(6):2878–2892. https://​doi.​org/​10.​1002/​jgrb.​
50219
Takeo A, Kawakatsu H, Isse T, Nishida K, Sugioka H, Ito A, Shiobara H,
Suetsugu D (2016) Seismic azimuthal anisotropy in the oceanic lithosphere and asthenosphere from broadband surface wave analysis of
OBS array records at 60 ma seafloor. J Geophys Res 121(3):1927–1947.
https://​doi.​org/​10.​1002/​2015J​B0124​29
Takeo A, Forsyth DW, Weeraratne DS, Nishida K (2014) Estimation of azimuthal anisotropy in the NW pacific from seismic ambient noise in
seafloor records. Geophys J Int 199(1):11–22. https://​doi.​org/​10.​1093/​
gji/​ggu240
Takeo A, Kawakatsu H, Isse T, Nishida K, Shiobara H, Sugioka H, Ito A, Utada
H (2018) In situ characterization of the Lithosphere–Asthenospheresystem beneath NW pacific ocean via broadband dispersion survey
with two OBS arrays. Geochem Geophys Geosyst. https://​doi.​org/​10.​1029/​2018G​C0075​88
Takeo A, Nishida K, Aoyama H, Ishise M, Kai T, Kurihara R, Maeda T, Mizutani Y, Nakashima Y, Nagahara S, Wang X, Ye L, Akuhara T, Aoki Y(2022) S-wave modelling of the Showa-Shinzan lava dome in Usu
Volcano, northern Japan, from seismic observations. Geophys J Int
230(3):1662–1678. https://​doi.​org/​10.​1093/​gji/​ggac1​11
Tanimoto T (1990) Modelling curved surface wave paths: membrane surface
wave synthetics. Geophys J Int 102:89–100. https://​doi.​org/​10.​1111/j.​
1365-​246X.​1990.​tb005​32.x
Tanimoto T, Anderson DL (1985) Lateral heterogeneity and azimuthal anisotropy of the upper mantle: Love and Rayleigh waves 100–250 s. J
Geophys Res 90(B2):1842. https://​doi.​org/​10.​1029/​jb090​ib02p​01842
Tanimoto T, Um J (1999) Cause of continuous oscillations of the earth. J
Geophys Res 104(B12):28723–28739. https://​doi.​org/​10.​1029/​1999J​
B9002​52
Tarantola A (2005) Inverse problem theory and methods for model parameter
estimation. Other titles in applied mathematics. Society for Industrial
and Applied Mathematics, Philadelphia. https://​doi.​org/​10.​1137/1.​
97808​98717​921
Tarantola A, Valette B (1982) Generalized nonlinear inverse problems solved
using the least squares criterion. Rev Geophys 20(2):219. https://​doi.​
org/​10.​1029/​RG020​i002p​00219
Tauzin B, Pham T-S, Tkalčić H (2019) Receiver functions from seismic interferometry: a practical guide. Geophys J Int 217(1):1–24. https://​doi.​org/​10.​1093/​gji/​ggz002
Thurber C, Ritsema J (2015) Theory and observations-seismic tomography andinverse methods. Elsevier, Amsterdam, pp 307–337. https://​doi.​org/​10.​1016/​B978-0-​444-​53802-4.​00009-9
Tian D, Uieda L, Leong WJ, Schlitzer W, FroÅNhlich Y, Grund M, Jones M, Toney
L, Yao J, Magen Y, Tong 1924 J-H, Materna K, Belem A, Newton T, AnantA, Ziebarth M, Quinn J, Wessel P (2023) PyGMT: A Python 1925 Interface
for the Generic Mapping Tools. https://​doi.​org/​10.​5281/​zenodo.​83031​
86
Tokimatsu K, Tamura S, Kojima H (1992) Effects of multiple modes on Rayleighwave dispersion characteristics. J Geotech Eng 118(10):1529–1543.
https://​doi.​org/​10.​1061/​(asce)​0733-​9410(1992)​118:​10(1529)
Tomar G, Shapiro NM, Mordret A, Singh SC, Montagner J-P (2017) Radialanisotropy in Valhall: ambient noise-based studies of Scholte and Lovewaves. Geophys J Int 208(3):1524–1539. https://​doi.​org/​10.​1093/​gji/​ggw480Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Tomar G, Stutzmann E, Mordret A, Montagner J-P, Singh SC, Shapiro NM (2018)
Joint inversion of the first overtone and fundamental mode for deep
imaging at the Valhall oil field using ambient noise. Geophys J Int
214(1):122–132. https://​doi.​org/​10.​1093/​gji/​ggy122
Trabant C, Hutko AR, Bahavar M, Karstens R, Ahern T, Aster R (2012) Data
products at the IRIS DMC: stepping stones for research and otherapplications. Seismol Res Lett 83(5):846–854. https://​doi.​org/​10.​1785/​
02201​20032
Tribaldos RV, Ajo-Franklin JB (2021) Aquifer monitoring using ambient seismic
noise recorded with distributed acoustic sensing (DAS) deployed ondark fiber. J Geophys Res. https://​doi.​org/​10.​1029/​2020j​b0210​04
Tromp J, Dahlen FA (1992) Variational principles for surface wave propagation
on a laterally heterogeneous Earth—II. Frequency-domain JWKB theory.
Geophys J Int 109(3):599–619. https://​doi.​org/​10.​1111/j.​1365-​246X.​
1992.​tb001​20.x
Tromp J, Fa Dahlen (1992) Variational principles for surface wave propagation on a laterally heterogeneous Earth-I. Time-domain JWKB theory.
Geophys J Int 109(3):581–598. https://​doi.​org/​10.​1111/j.​1365-​246X.​
1992.​tb001​19.x
Tromp J, Fa Dahlen (1993) Variational principles for surface wave propagation on a laterally heterogeneous Earth-III. Potential representation.
Geophys J Int 112(2):195–209. https://​doi.​org/​10.​1111/j.​1365-​246X.​
1993.​tb014​49.x
Tromp J, Luo Y, Hanasoge S, Peter D (2010) Noise cross-correlation sensitivity
kernels. Geophys J Int 183(2):791–819. https://​doi.​org/​10.​1111/j.​1365-​
246X.​2010.​04721.x
Tsai VC (2011) Understanding the amplitudes of noise correlation measurements. J Geophys Res 116(B9):09311. https://​doi.​org/​10.​1029/​2011J​
B0084​83
Tsai VC, Moschetti MP (2010) An explicit relationship between time-domainnoise correlation and spatial autocorrelation (SPAC) results. Geophys JInt. https://​doi.​org/​10.​1111/j.​1365-​246X.​2010.​04633.x
Ueno T, Saito T, Shiomi K, Haryu Y (2015) Monitoring the instrument responseof the high-sensitivity seismograph network in Japan (Hi-net): effects
of response changes on seismic interferometry analysis. Earth Planets
Space 67(1):135. https://​doi.​org/​10.​1186/​s40623-​015-​0305-0
van der Walt S, Colbert SC, Varoquaux G (2011) The NumPy array: a structure
for efficient numerical computation. Comput Sci Eng 13(2):22–30.
https://​doi.​org/​10.​1109/​mcse.​2011.​37
Viens L, Perton M, Spica ZJ, Nishida K, Yamada T, Shinohara M (2022) Understanding surface wave modal content for high-resolution imaging ofsubmarine sediments with distributed acoustic sensing. Geophys J Int
232(3):1668–1683. https://​doi.​org/​10.​1093/​gji/​ggac4​20
Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D,
Burovski E, Peterson P, Weckesser W, Bright J, van der Walt SJ, Brett M,
Wilson J, Millman KJ, Mayorov N, Nelson ARJ, Jones E, Kern R, Larson E,
Carey CJ, Polat I, Feng Y, Moore EW, VanderPlas J, Laxalde D, Perktold J,
Cimrman R, Henriksen I, Quintero EA, Harris CR, Archibald AM, Ribeiro
AH, Pedregosa F, van Mulbregt P, SciPy 1.0 Contributors (2020) SciPy 1.0:
fundamental algorithms for scientific computing in Python. Nat Methods 17(3):261–272. https://​doi.​org/​10.​1038/​s41592-​019-​0686-2
Wang J, Wu G, Chen X (2019) Frequency-Bessel transform method for effective
imaging of higher-mode Rayleigh dispersion curves from ambient
seismic noise data. J Geophys Res 124(4):3708–3723. https://​doi.​org/​10.​1029/​2018j​b0165​95
Wang K, Liu Q, Yang Y (2019) Three-dimensional sensitivity kernels for
multicomponent empirical Green’s functions from ambient noise:
Methodology and application to adjoint tomography. J Geophys Res
124(6):5794–5810. https://​doi.​org/​10.​1029/​2018j​b0170​20
Wang Y, Lin F-C, Schmandt B, Farrell J (2017) Ambient noise tomography
across mount St. Helens using a dense seismic array. J Geophys Res
122(6):4492–4508. https://​doi.​org/​10.​1002/​2016J​B0137​69
Wang K, Jiang C, Yang Y, Schulte-Pelkum V, Liu Q (2020) Crustal deformation
in southern California constrained by radial anisotropy from ambient
noise adjoint tomography. Geophys Res Lett. https://​doi.​org/​10.​1029/​
2020g​l0885​80
Wapenaar K, Draganov D, Snieder R, Campman X, Verdel A (2010) Tutorial
on seismic interferometry: part 1-basic principles and applications.
Geophysics 75(5):75–195

Page 41 of 42

Wapenaar K, Slob E, Snieder R, Curtis A (2010) Tutorial on seismic interferometry: part 2-underlying theory and new advances. Geophysics
75(5):75–211
Weaver RL (2011) On the amplitudes of correlations and the inference ofattenuations, specific intensities and site factors from ambient noise. C
R Geosci 343(8):615–622. https://​doi.​org/​10.​1016/j.​crte.​2011.​07.​001
Weaver R, Froment B, Campillo M (2009) On the correlation of non-isotropically
distributed ballistic scalar diffuse waves. J Acoust Soc Am 126(4):1817.
https://​doi.​org/​10.​1121/1.​32033​59
Webb SC (1998) Broadband seismology and noise under the ocean. Rev
Geophys 36(1):105–142
Wegler U, Sens-Schönfelder C, Sens-Schönfelder C (2007) Fault zone monitoring with passive image interferometry. Geophys J Int 168(3):1029–1033.
https://​doi.​org/​10.​1111/j.​1365-​246X.​2006.​03284.x
Wessel P, Luis JF, Uieda L, Scharroo R, Wobbe F, Smith WHF, Tian D (2019)
The Generic Mapping Tools Version 6. Geochem Geophys Geosyst
20(11):5556–64. https://​doi.​org/​10.​1029/​2019g​c0085​15
White MCA, Fang H, Nakata N, Ben-Zion Y (2020) PyKonal: a python package
for solving the Eikonal equation in spherical and cartesian coordinates
using the fast marching method. Seismol Res Lett. https://​doi.​org/​10.​
1785/​02201​90318
Wiener N (1947) Cybernetics: or control and communication in the animal andthe machine. MIT Press, Cambridge
Williams EF, Fernández-Ruiz MR, Magalhaes R, Vanthillo R, Zhan Z, GonzálezHerráez M, Martins HF (2021) Scholte wave inversion and passive
source imaging with ocean-bottom DAS. Lead Edge 40(8):576–583.
https://​doi.​org/​10.​1190/​tle40​080576.1
Woodhouse JH, Dziewonski AM (1984) Mapping the upper mantle: Threedimensional modeling of earth structure by inversion of seismic
waveforms. J Geophys Res 89(B7):5953–5986. https://​doi.​org/​10.​1029/​
jb089​ib07p​05953
Yamaya L, Mochizuki K, Akuhara T, Nishida K (2021) Sedimentary structure
derived from multi-mode ambient noise tomography with dense OBS
network at the Japan trench. J Geophys Res. https://​doi.​org/​10.​1029/​
2021j​b0217​89
Yang X, Gao H (2018) Full-wave seismic tomography in the northeastern
United States: new insights into the uplift mechanism of the Adirondack mountains. Geophys Res Lett 45(12):5992–6000. https://​doi.​org/​
10.​1029/​2018g​l0784​38
Yang Y, Ritzwoller MH, Levshin AL, Shapiro NM (2007) Ambient noise Rayleighwave tomography across Europe. Geophys J Int 168(1):259–274. https://​doi.​org/​10.​1111/j.​1365-​246X.​2006.​03203.x
Yang Y, Ritzwoller MH, Lin F-C, Moschetti MP, Shapiro NM (2008) Structure of
the crust and uppermost mantle beneath the western United States
revealed by ambient noise and earthquake tomography. J Geophys
Res. https://​doi.​org/​10.​1029/​2008j​b0058​33
Yang Y, Langston CA, Powell CA, Thomas WA (2022) Full waveform ambient
noise tomography for the northern Mississippi embayment. J Geophys
Res. https://​doi.​org/​10.​1029/​2021j​b0222​67
Yao H, van der Hilst RD, Montagner J-P (2010) Heterogeneity and anisotropyof the lithosphere of SE Tibet from surface wave array tomography. J
Geophys Res. https://​doi.​org/​10.​1029/​2009J​B0071​42
Yao H, Gouédard P, Collins JA, McGuire JJ, van der Hilst RD (2011) Structureof young east pacific rise lithosphere from ambient noise correlation
analysis of fundamental- and higher-mode Scholte-Rayleigh waves. C R
Geosci 343(8–9):571–583. https://​doi.​org/​10.​1016/j.​crte.​2011.​04.​004
Yoshizawa K, Kennett BLN (2002) Non-linear waveform inversion for surface
waves with a neighbourhood algorithm-application to multimodedispersion measurements. Geophys J Int 149(1):118–133. https://​doi.​
org/​10.​1046/j.​1365-​246X.​2002.​01634.x
Yoshizawa K, Kennett BLN (2004) Multimode surface wave tomography for
the Australian region using a three-stage approach incorporating finite
frequency effects. J Geophys Res 109(B2):1–19. https://​doi.​org/​10.​1029/​
2002J​B0022​54
Yoshizawa K, Kennett BLN (2005) Sensitivity kernels for finite-frequency surface
waves. Geophys J Int 162(3):910–926. https://​doi.​org/​10.​1111/j.​1365-​
246X.​2005.​02707.x
Zagoskin A (2014) Quantum theory of many-body systems. Springer, New
York. https://​doi.​org/​10.​1007/​978-3-​319-​07049-0

Nishida et al. Progress in Earth and Planetary Science(2024) 11:4

Zhao K, Yang Y, Luo Y (2020) Broadband finite frequency ambient noise
tomography: a case study in the western United States using USArraystations. J Geophys Res. https://​doi.​org/​10.​1029/​2019j​b0193​14
Zheng S, Sun X, Song X, Yang Y, Ritzwoller MH (2008) Surface wave tomography of China from ambient seismic noise correlation. Geochem
Geophys Geosyst 9(5):1–8. https://​doi.​org/​10.​1029/​2008G​C0019​81
Zhou L, Song X, Weaver RL (2020) Retrieval of amplitude and attenuation fromambient seismic noise: synthetic data and practical considerations.
Geophys J Int 222(1):544–559. https://​doi.​org/​10.​1093/​gji/​ggaa1​94

Page 42 of 42