Supporting information for “Vertical slice
ocean tomography with seismic waves”
Jörn Callies
1
, Wenbo Wu
1
,
2
, Shirui Peng
1
, Zhongwen Zhan
1
1
California Institute of Technology, Pasadena, CA, USA
2
Wood Hole Oceanographic Institution, Woods Hole, MA, USA
Contents of this file
1. Texts S1 to S2
2. Figures S1 to S9
3. Table S1
Corresponding author: Jörn Callies, jcallies@caltech.edu
1
Text S1 Inversion for travel time anomalies
The procedure to infer travel time anomalies relative to an unknown but common reference
from measurements of
T
-wave arrival time changes and origin time corrections (relative
to a cataloged event time) is similar to that described in Wu et al. (2020), but here we
improve on that inversion in a number of ways. In particular, we invert for the
T
-wave delays
5
and origin-time corrections separately, which allows us to more naturally use multiple
T
-
wave and land stations as well as a set of
T
-wave frequencies. Instead of the somewhat
ad hoc
regularization employed in Wu et al. (2020), we here impose a full set of prior and
noise statistics and use a Gauss–Markov estimator. See, for example, Wunsch (2006) for an
introduction to the methods employed here.
10
Our primary goal is to infer time series of the travel time anomalies at the observed
frequencies
ω
1
, . . . ,
ω
l
, stacked into the vector
τ
=
τ
1
.
.
.
τ
l
.
(1)
These travel time anomalies, ascribed to changes in the ocean’s temperature field, are dif-
ferences between
T
-wave arrival anomalies and origin time corrections:
τ
=
D
a
with
D
=
I
m
−
I
m
.
.
.
.
.
.
I
m
−
I
m
and
a
=
a
(
T
)
1
.
.
.
a
(
T
)
l
a
(
L
)
,
(2)
where
a
(
T
)
i
contains the
T
-wave travel time anomalies for each event time at the different
15
frequencies, and
a
(
L
)
contains the origin time corrections inferred from the land stations.
Here,
I
m
denotes the identity matrix of size
m
, the number of unique events. All delays are
referenced to predicted arrival times based on the cataloged event times. The predictions
are made using the PREM solid earth model for the land stations and a constant
T
-wave
reference velocity of 1
.
5 km s
−
1
.
20
The travel time anomalies
τ
are related to the range-averaged temperature anomalies
through the sensitivity kernel:
τ
=
(
K
⊗
I
m
)
T
∆
z
,
where
K
=
K
1
T
.
.
.
K
l
T
and
T
=
T
1
.
.
.
T
n
(3)
are the matrices containing the range-integrated kernels as rows and a stack of the time
series of range-averaged temperature anomalies at the
n
vertical levels, respectively, and
⊗
denotes a Kronecker product. As described in the main text, we perform an SVD of the sen-
25
sitivity matrix
K
=
U
Λ
V
T
. We normalize such that
U
T
U
=
I
3
and
h
−
1
V
T
V
∆
z
=
I
l
, where
2
∆
z
=
50 m is the vertical grid spacing and
h
=
5 km is a reference depth. With this normal-
ization, the singular vectors are unitless and do not depend on the discretization (as long as
it is fine enough). We denote the projections of the temperature anomalies onto the singular
vectors by
c
=
h
−
1
(
V
T
⊗
I
m
)
T
∆
z
, so that
τ
=
h
(
U
Λ
⊗
I
m
)
c
, which can be inverted for the time
30
series
c
=
h
−
1
(
Λ
−
1
U
T
⊗
I
m
)
τ
.
We cannot, however, observe the travel time anomalies directly. We only observe
changes
in the travel time between repeating earthquakes. The design matrix thus takes differences
of travel time anomalies between event pairs:
E
=
μ
I
l
⊗
X
(
T
)
0
0
X
(
L
)
¶
.
(4)
The pair matrices
X
(
T
)
and
X
(
L
)
take these differences between events at all
T
-wave and
35
land stations involved. For example, if there were four pairs connecting five events observed
at one
T
-wave section, we might have
X
(
T
)
=
−
1
1
−
1
1
−
1
1
−
1
1
.
(5)
If pairs 1, 2, and 4 are observed at land station 1, and pairs 1 and 3 are observed on land
station 2, the land station pair matrix would be
X
(
L
)
=
−
1
1
−
1
1
−
1
1
−
1
1
−
1
1
.
(6)
The observed arrival time offsets, again measured relative to the cataloged event times, are
40
then
δ
=
E
a
+
n
with
δ
=
δ
(
T
)
1
.
.
.
δ
(
T
)
l
δ
(
L
)
(7)
containing a stack of the
T
-wave offsets measured at different frequencies and the land
offsets. The observational noise is denoted by
n
.
We specify prior statistics for the projections onto the singular vectors (i.e., for
c
), and
we assign these covariances to the
T
-wave anomalies
a
(
T
)
1
, . . . ,
a
(
T
)
l
:
45
R
=
μ
h
2
(
U
Λ
⊗
I
m
)(
Σ
⊗
C
)(
U
Λ
⊗
I
m
)
T
0
0
0
¶
+
σ
2
1
l
+
1
⊗
I
m
,
(8)
3
where
Σ
is a diagonal matrix that specified the prior variances in the time series of the
projections onto the singular vectors (Table S1). The diagonal nature of
Σ
means that we
assign no prior correlations between the projections onto the singular vectors. The matrix
C
specifies the time correlation
C
i j
=
e
−|
t
i
−
t
j
|
/
λ
, which we assume to be identical for all singular
vectors. The expected standard deviation of origin time corrections is specified as
σ
=
5 s, and
50
1
l
+
1
denotes a square matrix of size
l
+
1 filled with ones. This prescribes perfect correlation
between the origin time correction and the part of the
T
-wave anomalies that is due to origin
time errors.
In addition to the random process with exponential correlation function, we also include
deterministic annual and semi-annual cycles as well as a linear trend. We do so by append-
55
ing their coefficients to the vector
a
and extending the matrices
E
,
R
, and
D
. For the design
matrix,
E
→
μ
A
E
0
¶
(9)
with
A
=
¡
I
l
⊗
X
(
T
)
t
I
l
⊗
X
(
T
)
cos
ω
t
I
l
⊗
X
(
T
)
sin
ω
t
I
l
⊗
X
(
T
)
cos 2
ω
t
I
l
⊗
X
(
T
)
sin 2
ω
t
¢
,
(10)
where
t
contains the event times referenced to the midpoint between the first event and last
event. For the prior covariance matrix,
60
R
→
μ
R
0
0
h
2
(
I
5
⊗
U
Λ
)
Ξ
(
I
5
⊗
U
Λ
)
T
¶
,
(11)
where
Ξ
is a diagonal matrix containing the prior variances for the trend or seasonal am-
plitudes for the projections onto the singular vectors. Like for the random components, we
assign no prior correlation between the trends and seasonal amplitudes of these projections.
For the difference matrix,
D
→
¡
D B
¢
,
(12)
with
65
B
=
¡
X
(
T
)
t
X
(
T
)
cos
ω
t
X
(
T
)
sin
ω
t
X
(
T
)
cos 2
ω
t
X
(
T
)
sin 2
ω
t
¢
,
(13)
such that the trends and seasonal variations are included in the travel time anomalies
τ
.
Once we specify the noise statistics
N
=
σ
n
2
I
with
σ
n
=
0
.
01 s, we now have all informa-
tion required for a Gauss–Markov estimate:
̃
a
=
PE
T
N
−
1
δ
with
P
=
¡
R
−
1
+
E
T
N
−
1
E
¢
−
1
.
(14)
The desired estimate of the travel time anomalies is then
̃
τ
=
D
̃
a
, and the posterior covari-
ance
P
can be propagated to this estimate as
DPD
T
. Similarly, estimates of the projections
70
onto singular vectors
̃
c
, of the trends, and seasonal amplitudes of the signal can be obtained
from
̃
a
, and their uncertainties can be calculated from
P
.
The prescription
σ
n
=
0
.
01 s is meant to capture errors due to the changes in the source
location between repeating events (Wu et al., 2020), due to noise affecting the cross-correlation
4
corr. time
random
trend
12 mo.
6 mo.
(days)
(mK)
(mK yr
−
1
)
(mK)
(mK)
Diego Garcia
15
15 10
5
4
2
1
15 10
5
15 10
5
Cape Leeuwin
30
20
5
5
4
2
1
20
5
5
5
2
2
Table S1: Parameters of the prior covariances of the projections onto the singular vectors.
Where three parameters are given, they are for the three singular vectors. For the seasonal
signals, the cos and sin terms are each assigned half the indicated prior variance.
function, and due to hydrophone motion. The moored CTBTO hydrophones are displaced by
75
local currents, which contributes measurement error because the moorings are not navi-
gated. Nichols and Bradley (2017) estimated that the CTBTO hydrophones at Wake Island,
similar in design to the ones used here, can be displaced by 1/20 of the length of their moor-
ing lines. The mooring lines for H08S2 and H01W2 are 562 m and 570 m long, respectively,
producing maximum horizontal displacements of less than 30 m. A horizontal displacement
80
of 30 m corresponds to a travel time anomaly of 0
.
02 s, which is equal to 2
σ
n
. A more re-
fined prescription of
N
might split the errors into their contributing factors and take the
respective error correlations into account. We leave this refinement to future work.
Text S2 Cycle-skipping correction
Acoustic modes are dispersive. If the sounds speed profile is perturbed, a dispersive mode’s
85
phase and group will shift by different amounts. If the differential shift between a mode’s
phase and group is large enough, the correlation function will peak not at the correct phase
shift but at the next peak over. The modal dispersion causes a cycle skip.
We here illustrate this process using a simple example of a Gaussian wave packet prop-
agating through Munk’s canonical sound speed profile. The observed
T
-wave delays show a
90
tell-tale sign of this process, suggesting that the same process occurs in the real system. The
simple physics causing the cycle skipping allows for a robust identification and correction
procedure, described below.
To illustrate the process causing cycle skipping, we consider modal propagation through
a range-independent ocean with slowness profile
S
(
z
), following Munk et al. (1995). A signal
95
propagating a distance
L
has a phase travel time
τ
=
Ls
, where
s
=
k
/
ω
is the modal phase
slowness that equals the slowness
S
(
z
) at the turning depths. The modal structure
P
is
defined by
d
2
P
d
z
2
+
¡
ω
2
S
2
−
k
2
¢
P
=
0
(15)
with boundary conditions
P
=
0 at
z
=
0 (the surface) and d
P
/d
z
=
0 at
z
=−
h
(the bottom).
We apply the normalization
R
d
z P
2
=
h
. The group slowness is
100
s
g
=
d
k
d
ω
=
1
sh
Z
d
z S
2
P
2
,
(16)
5
where the integration is over the full depth. The mode’s group travel time is
τ
g
=
Ls
g
.
We now perturb the slowness profile
S
(
z
) to
S
(
z
)
+
∆
S
(
z
). The difference in the modal
phase slowness between the reference and perturbed profiles is
∆
s
=
1
sh
Z
d
z P
2
S
∆
S
.
(17)
The change in group slowness, in contrast, is
∆
s
g
=
1
sh
Z
d
z
h³
2
−
s
g
s
́
P
2
+
ω
(
P
2
)
ω
i
S
∆
S
,
(18)
where (
P
2
)
ω
=
2
PP
ω
, and
P
ω
is defined by
105
d
2
P
ω
d
z
2
+
¡
ω
2
S
2
−
k
2
¢
P
ω
=−
2
¡
ω
S
2
−
ks
g
¢
P
(19)
with the same boundary conditions as for
P
and the additional constraint that
R
d
z P
ω
=
0.
Because variations in the slowness are much smaller than their mean values, the sensitivity
profiles are roughly
P
2
/
h
for the phase and [
P
2
+
ω
(
P
2
)
ω
]/
h
for the group.
These phase and group sensitivity profiles thus have different structures (Fig. S2). We
illustrate their behavior using the lowest mode at
ω
/2
π
=
2
.
5 Hz for the canonical temperate
110
profile of Munk (1974) and Munk et al. (1995; eq. 2.2.13). Above 1
.
6 km depth, the group
sensitivity to sound speed perturbations is larger than the phase sensitivity. If the temper-
ature anomalies that cause these sound speed perturbations are dominated by this depth
range, as is typical, the group delay
∆
τ
g
=
L
∆
s
g
will exceed the phase delay
∆
τ
=
L
∆
s
. If the
group delay is more than a half period larger than the phase delay, cycle skipping is likely.
115
In this scenario, cycle skipping is always forward: the skip goes in the same direction as the
delay. It is important to correct this cycle skipping because otherwise large anomalies would
be overestimated systematically.
As an example, consider a slowness perturbation that is caused by a single low dynamical
mode (Fig. S2). For the buoyancy frequency profile
N
=
N
0
e
z
/
d
underlying the canonical
120
sound speed profile, the sound speed perturbation is proportional to
N
2
G
. The dynamical
mode
G
for vertical displacement is defined by
d
2
G
d
z
2
+
N
2
c
2
G
=
0
(20)
with
G
=
0 at
z
=
0 and
z
= −
h
, and a normalization
R
d
z N
2
G
2
=
h
is applied. We use
h
=
5 km, but this has little impact on the results because neither the first acoustic mode
nor the first dynamical mode has much amplitude in the bottom 1 km or so. For the first
125
dynamical,
∆
τ
g
/
∆
τ
=
1
.
61, so cycle skipping is likely for
|
∆
τ
|>
0
.
33 s. This is about where we
observe cycle skipping to occur in the real
T
waveforms received at Cape Leeuwin.
To further illustrate the cycle skipping, we consider a signal with the source function
e
−
σ
2
t
2
cos
ω
t
. We again use the frequency
ω
/2
π
=
2
.
5 Hz, and we impose a bandwidth
σ
/2
π
=
0
.
5 Hz corresponding to the filtering we use for the real data. We then propagate this signal
130
6
0
1
2
3
4
5
depth (km)
Diego Garcia
(a)
SEM
2.00 Hz
3.00 Hz
4.00 Hz
Diego Garcia
(b)
modal
6
4
2
0
sensitivity (s K
1
km
1
)
0
1
2
3
4
5
depth (km)
Cape Leeuwin
(c)
SEM
2.50 Hz
3.25 Hz
4.00 Hz
6
4
2
0
sensitivity (s K
1
km
1
)
Cape Leeuwin
(d)
modal
Figure S1: Sensitivity kernels calculated using two-dimensional numerical calculations as
well as assuming the propagation of the first acoustic mode through a range-independent
ocean. These kernels are shown for all considered frequencies and for both paths.
assuming a dispersion relation linearized around
ω
. We impose phase lags
∆
τ
=
0
.
0 s
,
0
.
3 s
,
and 0
.
6 s and group lags
∆
τ
g
=
1
.
61
∆
τ
, corresponding to propagation through an anomaly
caused by the first dynamical mode (Fig. S3). Correlating the delayed signals with the signal
without lag, the signal with
∆
τ
=
0
.
3 s produces a maximum of the correlation function at the
correct phase lag. But the signal with
∆
τ
=
0
.
6 s experiences cycle skipping: the maximum
135
of the correlation function shifts by one period to a lag of 1
.
0 s. The secondary maximum to
the left is at the correct lag. If cycle skipping can be identified, it can easily be corrected by
picking the adjacent peak in the correlation function.
This type of cycle skipping caused by dispersive modal propagation can be identified in
the real data. If the first acoustic mode dominates and temperature anomalies are surface-
140
intensified, the cycle skipping causes an increase in the magnitude of the low-frequency (say,
2
.
5 Hz) phase delay. At the same time, it causes a decrease in the differential delay between
a higher frequency (say, 4 Hz) and this low frequency. Plotting the measured low-frequency
and differential delays against one another reveals three distinct clusters: a center cluster
with small low-frequency delays that is not affected by cycle skipping and adjacent clus-
145
ters with larger-magnitude low-frequency delays that suffer from forward or backward cycle
skipping (Fig. S4, S5).
7
0.000
0.002
0.004
buoyancy freq. (s
1
)
0
1
2
3
4
5
depth (km)
stratication
(a)
0.65
0.66
slowness (s km
1
)
sound slowness
(b)
2
0
2
mode
P
acoustic modes
(c)
0.0
2.5
5.0
sensitivity (km
1
)
sensitivity
(d)
phase
group
0.00
0.01
mode
N
2
G
(s
2
)
dyn. modes
(e)
Figure S2: Illustration of dispersive propagation of low-frequency acoustic modes. Shown
are profiles of (a) the buoyancy frequency
N
, (b) Munk’s canonical slowness profile
S
(blue)
and the phase slownesses
s
of the first acoustic mode (solid orange) and two higher modes
(dotted orange), (c) the first acoustic mode (solid) and two higher modes (dotted), (d) the
phase and group sensitivities
P
2
S
/
sh
(green) and [(2
−
s
g
/
s
)
P
2
+
ω
(
P
2
)
ω
]
S
/
sh
(red) for the
first acoustic mode, and (e) the first dynamical mode (solid) and two higher modes (dotted).
We identify pairs in these clusters by employing a Gaussian mixture model with three
components that have shared covariances. We correct pairs in the left cluster by shifting
to the next maximum of the correlation function to the right of the original maximum, and
150
we correct pairs in the right cluster by shifting to the next maximum to the left. We then
probe for additional cycle skips (or reversals of corrections applied after the cluster analysis)
by looking for reductions in the cost of the inversion for travel time anomalies (cf. Wu et
al., 2020). We allow such additional corrections if the Gaussian mixture model indicates a
probability greater than 0.1 % to belong to a cluster different from that identified as most
155
likely. For the path to Diego Garcia, this procedure corrects 40 pairs to the left and 56 pairs
to the right (out of a total of 3831 used pairs). For the path to Cape Leeuwin, cycle skipping
is more common, both because the path is longer and because travel time anomalies are
somewhat larger. Here, 230 pairs are corrected to the left, and 393 pairs are corrected to the
right (out of a total of 3032 used pairs).
160
References
Munk, W., P. Worcester, and C. Wunsch (1995)
Ocean Acoustic Tomography
. Cambridge Uni-
versity Press, p. 433.
Munk, W. H. (1974) Sound channel in an exponentially stratified ocean, with application to
SOFAR.
Journal of the Acoustical Society of America
55 (2), 220–226.
DOI
: 10.1121/1.
165
1914492.
8
1.0
0.5
0.0
0.5
1.0
signal
= 0.0 s
= 0.3 s
= 0.6 s
2
1
0
1
2
3
lag (s)
1.0
0.5
0.0
0.5
1.0
cross-correlation
Figure S3: Illustration of cycle skipping using synthetic signals. The top panel shows
three signals that have experienced different amounts of phase lags
∆
τ
and group lags
∆
τ
g
=
1
.
61
∆
τ
. The triangles trace the peak that is at the center of the signal without
lag, illustrating the phase shift. The bottom panel shows the cross-correlation functions
between the lagged signals and the reference signal without lag. The maximum of the cross-
correlation function occurs at the correct phase lag (marked by triangle) for the signal with
∆
τ
=
0
.
3 s. Cycle skipping occurs for the signal with
∆
τ
=
0
.
6 s: the maximum (marked by
circle) is offset by one period from the correct lag, which coincides with a secondary maxi-
mum (marked by triangle)
Nichols, S. M. and D. L. Bradley (Oct. 2017) In Situ Shape Estimation of Triangular Moored
Hydrophone Arrays Using Ambient Signals.
IEEE Journal of Oceanic Engineering
42 (4),
923–935.
DOI
: 10.1109/JOE.2016.2625378.
Wu, W., Z. Zhan, S. Peng, S. Ni, and J. Callies (2020) Seismic ocean thermometry.
Science
170
1515 (6510), 1510–1515.
DOI
: 10.1126/science.abb9519.
Wunsch, C. (2006)
Discrete Inverse and State Estimation Problems
. Cambridge University
Press, p. 371.
9
1.0
0.5
0.0
0.5
1.0
low-frequency delay (s)
0.20
0.15
0.10
0.05
0.00
0.05
0.10
0.15
0.20
differential delay (s)
after clustering
(a)
1.0
0.5
0.0
0.5
1.0
low-frequency delay (s)
after additional corrections
(b)
Figure S4: Cycle-skipping corrections for the path to Diego Garcia. Shown are the low-
frequency delays (2
.
0 Hz) vs. the differential delays (4
.
0 Hz minus 2
.
0 Hz). The green dots
indicate identified cycle skipping that is corrected to the right, the blue dots indicate cycle
skipping corrected to the left, and orange dots indicate pairs not identified as being affected
by cycle skipping. The identifications are (a) after the cluster analysis and (b) after addition
corrections (or reversals) based on the inversion cost.
10