Supporting information for “Seismic ocean
thermometry using CTBTO hydrophones”
Wenbo Wu
1,2
, Zhichao Shen
1,2
, Shirui Peng
2
, Zhongwen Zhan
2
, and Jörn
Callies
2
1
Woods Hole Oceanographic Institution, Woods Hole, MA, USA
2
California Institute of Technology, Pasadena, CA, USA
Contents of this file
1. Text S1 to S3
2. Figures S1 to S8
Corresponding author: Wenbo Wu, wenbo.wu@whoi.edu
1
Text S1 Finding missing earthquakes with template matching
In this study, we apply the Template Matching (TM) method (see e.g., Anstey, 1964; Harris,
1991) to the two reference stations (GSI and PSI) to search for the small earthquakes miss-
ing in the current ISC catalog. The station PSI has been used in our previous study (Wu
et al., 2020). This station has a relatively long history starting before the 2005 M 8.6 earth-
5
quake, low noise levels and appropriate distances to the sources, all of which are desirable
for SOT purposes. We add another station GSI, which is the only publicly available station
within the source region and therefore very sensitive to small earthquakes. We download
the continuous waveforms of these two stations, remove the instrument responses and filter
them to 1.0–3.0 Hz. The 5374 ISC earthquakes in 2000–2016 around Nias Island (Fig. 1a in
10
the main text) are used to produce templates. For each ISC earthquake, its template seis-
mogram at a reference station is created by cutting a time window from 3 s before to 47 s
after the
P
-wave arrival. The 1-D Earth model PREM (Dziewonski and Anderson, 1981) is
used for the
P
-wave arrival prediction. Then the template matching is conducted by carry-
ing out waveform cross-correlation (CC) between the template seismogram and continuous
15
data in a moving window. A new earthquake is identified if the CC coefficient surpasses 0.6.
Scanning the continuous data allows us to identify all earthquakes that produce sufficiently
similar waveforms as the template. This process is repeated for each ISC earthquake and
each reference station.
TM produces a new catalog for each combination of ISC events and reference stations.
20
We check each of these catalogs and discard some obvious false detections. For example,
some templates have bad waveforms due to instrument glitches, and we opt to not use these
catalogs. Earthquakes above M5 are also discarded, since it is unlikely that such a big event
is missing in the ISC catalog. Then we merge the catalogs by eliminating redundant events
and thus obtain what we refer to as the “TM catalog”. Our elimination of redundant events
25
is simply based on their origin times—any two detections with an origin time difference
smaller than 6 s is identified as the same event. The origin time and location of the events
to be merged are assigned the same as the associated template having the maximum CC co-
efficient. The magnitude of a new detected event is estimated using its waveform amplitude
relative to the template.
30
This procedure expands the original catalog of 5374 ISC earthquakes to about 500 000
events (Fig. 2 in the main text). We note that this TM catalog is still incomplete, and there
must be many earthquakes that remain unidentified. Missing an earthquake could be due
to poor-quality data (i.e., low SNR or data gaps at the reference stations), or a lack of a good
template (i.e., no any ISC event has sufficiently similar waveform as the missing earth-
35
quake). The incompleteness is indicated by the clear earthquake shortage—for example,
the number of M
<
4 earthquakes are well below the Gutenberg–Richter curve. Thus, there
is still a potential to find more repeaters. We also note that some difficultly identified false
detections likely remain in the current catalog. For example, the number of earthquakes
with M4.5–5.0 are well above the Gutenberg–Richter prediction, which may be due to false
40
detections. Such false detections will not affect our SOT results, however, because they can-
not produce
T
-waves that pass our CC criteria at the hydrophone or
T
-wave station.
2
Text S2 Finding repeaters and measuring travel time changes
We search for repeaters by conducting waveform CC of both
P
and
T
waves between paired
earthquakes from the TM catalog. The
T
-wave arrival time change and relative origin times
45
between repeaters are simultaneously measured during the CC computations. We note that
one can identify the repeaters of an ISC event based on the CC results from the previous
step of CC between the
P
wave of a ISC event and the continuous waveforms. New repeating
pairs, however, could form between two non-ISC events. Thus, we redo the CC computations
using the TM catalog to obtain a complete list of repeating pairs.
50
For
P
waves at the four reference stations GSI, PSI, KUM, and WRAB, a 1.5–2.5 Hz or
1.0–3.0 Hz Butterworth band-pass filter is applied to the waveforms to improve the SNR.
Then we cut the
P
wave using a 50 s long time window for the CC computations. The time
window starts at 3 s before the first
P
-wave arrival, and this arrival time is predicted based
on PREM (Dziewonski and Anderson, 1981). In order to reduce computational cost, we only
55
search earthquake pairs whose cataloged locations are at most 60 km apart. A maximum
time shift up to 8 s is allowed in seeking the highest CC coefficient between paired earth-
quakes. Mathematically, the
P
-wave offset measurement may be written as
δ
L
=
arg max
δ
t
xcorr
¡
s
A
(
t
,
O
A
)
w
(
t
,
T
PREM
)
,
s
B
(
t
+
δ
t
,
O
B
)
w
(
t
,
T
PREM
)
¢
(1)
where xcorr denotes the normalized cross correlation between two windowed
P
-wave seis-
mograms
s
A
from event
A
and
s
B
from the paired event
B
. The two events’ cataloged origin
60
times are represented as
O
A
and
O
B
, respectively. The 50 s long window function
w
(
t
) starts
at 3 s before the predicted
P
-wave arrival. The PREM
P
-wave arrival prediction
T
P RE M
is
calculated using the cataloged source locations.
The
T
-wave CC calculation is similar as for
P
waves, but some parameters are different.
A frequency-domain Gaussian filter with a width of 0.5 Hz, rather than a Butterworth band-
65
pass filter, is applied to the
T
wave. In the accompanying paper (Callies et al., 2023), we
show that a
T
-wave travel time change is expected to depend on frequency. Compared to a
Butterworth filter, the Gaussian filter helps produce a smooth frequency dependence of the
measurements while concentrating the measurement on the target frequency. We choose a
center frequency of the Gaussian filter as low as possible while not too low to suffer from
70
high noise—2.5 Hz for H01 and 2.0 Hz for H08 and DGAR worked well in practice. The
path to H01 off Cape Leeuwin encounters shallow bathymetric features blocking the low-
frequency
T
wave, so it is not surprising that the 2.0 Hz data is noisier there and a higher
frequency is needed to increase the SNR. We measure the
T
-wave travel time change using a
60 s long time window that starts 10 s after the predicted
T
-wave arrival time. The
T
-wave
75
arrival time is calculated using a constant sound speed of 1
.
51 km s
−
1
and the cataloged
source locations. Similar to equation 1,
T
-wave offset can be expressed as
δ
T
=
arg max
δ
t
xcorr
¡
s
A
(
t
,
O
A
)
w
(
t
,
T
T
wave
)
,
s
B
(
t
+
δ
t
,
O
B
)
w
(
t
,
T
T
wave
)
¢
.
(2)
Note that the CTBTO deployed six hydrophones off the Diego Garcia atoll—a group of
three hydrophones to the south and another three to the north. Each group forms a triangle-
3
shape network with a station distance of about 2 km. Compared to the southern tri-net, the
80
northern tri-net records
T
waves from Sumatra with a much lower SNR because of blocking
by the atoll. The three hydrophones of the southern tri-net record similar
T
waves, and we
choose H08S2 for this study. For Cape Leeuwin, we use the data from H01W3.
The first criterion for a repeating pair to be used is that the
P
-wave CC coefficient from
at least one reference station exceeds 0.9. Many of the pairs passing this criterion, however,
85
fail to meet the second one—that the
T
-wave CC coefficient exceeds 0.6. This will occur in
the following three situations, (1) the
T
-wave station has a data gap, (2) a low
T
-wave SNR
causes low
T
-wave CC, (3) the SNR is high, but the
T
waveforms between the events is not
similar. As mentioned in Text S1 above, false detections in the TM step (e.g., due to recurring
instrument glitches) have high
P
-wave similarities at the reference stations, but they will
90
not meet the
T
-wave CC criterion.
Using this procedure, we detect 4704 repeating pairs for H08, 2956 pairs for DGAR, and
3766 pairs for H01. The
P
- and
T
-wave arrival time changes of these pairs are incorporated
into a linear system of equation for each
T
-wave station, and time series are estimated using
a Gauss–Markov inversion (see Callies et al., 2023).
95
Text S3 Sensitivity kernel calculations
Text S3.1 SPECFEM2D numerical simulations
Wu et al. (2020) used the Spectral Element Method (SEM) software package SPECFEM2D
(Komatitsch and Tromp, 1999) to compute the
T
-wave sensitivity kernels. The simulations
in this study are to a large degree similar. For example, the structure model is composed of
100
an ocean layer and an underlying solid earth; the sound speed is calculated in the same
way using ECCO and the Gibbs Seawater Toolbox; the crust has background structure
with a
P
-wave speed
v
p
=
5800 m s
−
1
, an
S
-wave speed
v
s
=
3200 m s
−
1
, and a density
ρ
=
2600 kg m
−
3
; this background structure is perturbed by small-scale heterogeneities; and
realistic bathymetry along the paths is used. More details about the model setup are given
105
in Wu et al. (2020).
We make some minor changes here to further improve the sensitivity kernel compu-
tations. First, we add a sediment layer below the seafloor. The East Indian Ocean basin
is widely covered with low-velocity sediment, which affects the sea-seabed interactions of
T
waves, especially in shallow water. To take these effects into account, we use the global
110
sediment thickness data from Straume et al. (2019) to build up a sediment layer atop the
crust. The seismic velocity and density of sediment as functions of depth are calculated us-
ing existing empirical results (Hamilton, 1976, 1979), while the extremely low
v
s
at the top
of the sediment is cut off at 1500 m s
−
1
for numerical stability. Second, we replace the Butter-
worth bandpass filter used in Wu et al. (2020) with the Gaussian filter, as described above.
115
Accordingly, the adjoint source in SEM simulation needs to be modified for consistency. For
the travel time measured with waveform CC, the adjoint source is given by
f
†
i
(
x
,
z
,
t
)
=
1
N
r
w
r
(
̄
t
−
t
)
∂
t
s
′
i
(
x
r
,
z
r
,
̄
t
−
t
)
δ
(
x
−
x
r
)
δ
(
z
−
z
r
)
(3)
4
and
N
r
=
Z
̄
t
0
w
r
(
̄
t
−
t
)
s
′
i
(
x
r
,
z
r
,
t
)
∂
tt
s
i
(
x
r
,
z
r
,
t
) d
t
,
(4)
where
s
i
(
x
r
,
z
r
,
t
) is the synthetic seismogram,
s
′
i
(
x
r
,
z
r
,
t
) is the filtered
s
i
(
x
r
,
z
r
,
t
),
w
r
(
t
) de-
notes the cross-correlation window, and
̄
t
is the seismogram length. The synthetic seismo-
120
gram is the vertical component of velocity for DGAR and pressure for H08 and H01. Because
we apply a Gaussian filter to the observed data, the adjoint source needs to be filtered in the
same way, i.e.,
s
′
i
(
x
r
,
z
r
,
ω
)
=
F
(
ω
)
F
∗
(
ω
)
s
(
ω
), where
F
(
ω
) is the Gaussian filter. Note that the
adjoint source needs to be filtered twice for zero phase shift to ensure consistency in timing
(Luo et al., 2013).
125
Text S3.2 Mode theory
For sound waves propagating in acoustic media, a travel time shift could be caused by small
perturbation to the media properties (i.e., density and sound speed). One may use the wave-
form cross-correlation scheme to measure the time lag of wave signals (i.e., pressure
p
as a
function of time
t
) in a given time window between
t
1
and
t
2
. This time lag should be equal
130
to (Dahlen et al., 2000)
∆
τ
=
R
t
2
t
1
∂
t
p
(
t
)
∆
p
(
t
) d
t
R
t
2
t
1
∂
tt
p
(
t
)
p
(
t
) d
t
,
(5)
where
∆
p
(
t
) is the pressure perturbation. In most applications, the data is filtered to some
frequency band between
ω
1
and
ω
2
before measuring
∆
τ
, so it is convenient to rewrite (5) in
the frequency domain as
∆
τ
=
Re
R
ω
2
ω
1
i
ω
ˆ
p
∗
(
ω
)
c
∆
p
(
ω
) d
ω
R
ω
2
ω
1
ω
2
|
ˆ
p
(
ω
)
|
2
d
ω
,
(6)
where carets denote Fourier transforms.
135
The acoustic perturbation in a waveguide that has a 1-D sound speed structure de-
pending on the depth
z
but not the range
x
can be expressed as a summation of range-
independent modes
p
(
ω
,
x
,
z
)
=
∞
X
m
=
1
A
m
(
ω
)
φ
m
(
ω
,
z
)
e
−
ik
m
x
,
(7)
where
m
is the mode number,
A
m
is amplitude of a complex number,
k
m
is the wavenumber,
and
φ
m
(
ω
,
z
) is the eigenfunction. The eigenfunctions satisfy the orthogonality condition
140
R
φ
m
(
z
)
φ
n
(
z
) d
z
=
δ
mn
. The
T
wave used in this study are very low frequency waves (i.e.,
<
4 Hz). In this frequency range, all modes higher than
m
=
1 have wide eigenfunctions
penetrating substantially into the solid Earth (Fig. S1), which causes complex interaction
with the seafloor. The consequence is that these higher modes are seriously attenuated,
leaving the fundamental mode (
m
=
1) as the dominant component of the
T
wave. Dropping
145
these low-energy higher modes results in a simplified version of (7):
p
(
ω
,
x
,
z
)
=
A
(
ω
)
φ
(
ω
,
z
)
e
−
ikx
,
(8)
5
where the order number
m
=
1 is omitted for convenience. With the assumption of small
sound speed perturbations and negligible effects from density changes, the first Born ap-
proximation is available for the pressure perturbations (e.g., Skarsoulis and Cornuelle,
2004; Skarsoulis et al., 2009). The signal change at a receiver located at (
x
r
,
z
r
) can be ob-
150
tained as
c
∆
p
(
ω
,
x
r
,
z
r
)
=
Ï
i
ω
2
∆
v
(
x
,
z
)
v
3
(
z
)
k
A
(
ω
)
φ
2
(
ω
,
z
)
φ
(
ω
,
z
r
)
e
−
ikx
r
d
x
d
z
.
(9)
Substituting (8) and (9) into (6), we obtain the travel time perturbation described by normal
modes. After some algebra, (6) becomes
∆
τ
=−
R
ω
2
ω
1
ω
2
|
A
(
ω
)
|
2
φ
2
(
ω
,
z
r
)
Î
∆
v
(
x
,
z
)
v
2
(
z
)
φ
2
(
ω
,
z
) d
x
d
z
d
ω
R
ω
2
ω
1
ω
2
|
A
(
ω
)
|
2
φ
2
(
ω
,
z
r
) d
ω
=
Ï
K
v
(
z
)
∆
v
(
x
,
z
)
v
(
z
)
d
x
d
z
,
(10)
with
K
v
(
z
)
=−
R
ω
2
ω
1
ω
2
|
A
(
ω
)
|
2
φ
2
(
ω
,
z
r
)
φ
2
(
ω
,
z
)
v
(
z
)
d
ω
R
ω
2
ω
1
ω
2
|
A
(
ω
)
|
2
φ
2
(
ω
,
z
r
) d
ω
.
(11)
The sensitivity kernel
K
v
(
z
) is a function of depth
z
but not range
x
, which is expected
155
given the assumption of a range-independent environment. We modify (11) by adding the
range dependence of the sound speed:
K
v
(
x
,
z
)
=−
R
ω
2
ω
1
ω
2
|
A
(
ω
)
|
2
φ
2
(
ω
,
x
r
,
z
r
)
φ
2
(
ω
,
x
,
z
)
v
(
x
,
z
)
d
ω
R
ω
2
ω
1
ω
2
|
A
(
ω
)
|
2
φ
2
(
ω
,
x
r
,
z
r
) d
ω
,
(12)
where the eigenfunctions can be computed using the local 1-D depth profile of sound speed.
Equation (12) allows us to calculate the sensitivity kernels for
T
waves arriving at DGAR,
H08, and H01. We extract the 2005–2015 mean temperature and salinity from the ECCO
160
state estimate (Forget et al., 2015) and calculate the sound speed using the Gibbs Seawa-
ter Toolbox (McDougall and Barker, 2011). The eigenfunctions
φ
are computed using the
Acoustic Tool software package (Porter, n.d.). The boundary conditions are specified with
pressure-free condition at the top and semi-infinite solid half space below the ocean layer
(Porter and Reiss, 1984, 1985). The solid media has a
P
-wave speed
v
p
=
4
.
0 km s
−
1
and a
165
density
ρ
=
3
.
0 g cm
−
1
. We assign Gaussian shapes to the amplitude
|
A
(
ω
)
|
, consistent with
the Gaussian filter applied to the observational data. The bathymetry is extracted from the
GEBCO data (http://www.gebco.net).
References
Anstey, N. A. (1964) Correlation techniques–a review.
Geophysical Prospecting
12 (4), 355–
170
382.
Callies, J., W. Wu, S. Peng, and Z. Zhan (2023) Vertical-Slice Ocean Tomography With Seis-
mic Waves.
Geophysical Research Letters
50 (8), e2023GL102881.
6
Dahlen, F., S.-H. Hung, and G. Nolet (2000) Fréchet kernels for finite-frequency travel-
times—I. Theory.
Geophysical Journal International
141 (1), 157–174.
175
Dziewonski, A. M. and D. L. Anderson (1981) Preliminary reference Earth model.
Physics of
the earth and planetary interiors
25 (4), 297–356.
Forget, G., J.-M. Campin, P. Heimbach, C. N. Hill, R. M. Ponte, and C. Wunsch (2015) ECCO
version 4: An integrated framework for non-linear inverse modeling and global ocean
state estimation.
Geoscientific Model Development
8 (10), 3071–3104.
180
Hamilton, E. L. (1976) Variations of density and porosity with depth in deep-sea sediments.
Journal of Sedimentary Research
46 (2), 280–300.
Hamilton, E. L. (1979) V p/V s and Poisson’s ratios in marine sediments and rocks.
The
Journal of the Acoustical Society of America
66 (4), 1093–1101.
Harris, D. B. (1991) A waveform correlation method for identifying quarry explosions.
Bul-
185
letin of the Seismological Society of America
81 (6), 2395–2418.
Komatitsch, D. and J. Tromp (1999) Introduction to the spectral element method for three-
dimensional seismic wave propagation.
Geophysical journal international
139 (3), 806–
822.
Luo, Y., J. Tromp, B. Denel, and H. Calandra (2013) 3D coupled acoustic-elastic migration
190
with topography and bathymetry based on spectral-element and adjoint methods.
Geo-
physics
78 (4), S193–S202.
McDougall, T. J. and P. M. Barker (2011) Getting started with TEOS-10 and the Gibbs Sea-
water (GSW) oceanographic toolbox.
Scor/Iapso WG
127, 1–28.
Porter, M. B. (n.d.) Acoustics Toolbox (). [Last viewed on 2-January-2023].
195
Porter, M. B. and E. L. Reiss (1984) A numerical method for ocean-acoustic normal modes.
The Journal of the Acoustical Society of America
76 (1), 244–252.
Porter, M. B. and E. L. Reiss (1985) A numerical method for bottom interacting ocean acous-
tic normal modes.
The Journal of the Acoustical Society of America
77 (5), 1760–1767.
Skarsoulis, E., B. Cornuelle, and M. Dzieciuch (2009) Travel-time sensitivity kernels in long-
200
range propagation.
The Journal of the Acoustical Society of America
126 (5), 2223–2233.
Skarsoulis, E. K. and B. D. Cornuelle (2004) Travel-time sensitivity kernels in ocean acoustic
tomography.
The Journal of the Acoustical Society of America
116 (1), 227–238.
Straume, E. O., C. Gaina, S. Medvedev, K. Hochmuth, K. Gohl, J. M. Whittaker, R. Ab-
dul Fattah, J. C. Doornenbal, and J. R. Hopper (2019) GlobSed: Updated total sediment
205
thickness in the world’s oceans.
Geochemistry, Geophysics, Geosystems
20 (4), 1756–1772.
Wu, W., Z. Zhan, S. Peng, S. Ni, and J. Callies (2020) Seismic ocean thermometry.
Science
369 (6510), 1510–1515.
7
Figure S1: The eigenfunctions of modes 1 and 2 for the Diego Garcia path. The corresponding
sound speed profile as a function of depth is computed using 2005–2017 mean temperature
and salinity from the ECCO state estimate (see Fig. 9b).
8
Figure S2: Sensitivity kernels
K
v
for 3.0 Hz
T
waves received at DGAR. (a)
K
v
computed
with SPECFEM2D. (b)
K
v
computed using mode theory. Only the fundamental mode is
included in the calculation. (c) Range-averaged
K
v
for the three segments 0–1000 km, 1000–
2000 km and 2000–2890 km.
9