Vertical-slice ocean tomography with seismic waves

Seismically generated sound waves that propagate through the ocean are used to infer temperature anomalies and their vertical structure in the deep East Indian Ocean. These T waves are generated by earthquakes off Sumatra and received by hydrophone stations off Diego Garcia and Cape Leeuwin. Between repeating earthquakes, a T wave's travel time changes in response to temperature anomalies along the wave's path. What part of the water column the travel time is sensitive to depends on the frequency of the wave, so measuring travel time changes at a few low frequencies constrains the vertical structure of the inferred temperature anomalies. These measurements reveal anomalies due to equatorial waves, mesoscale eddies, and decadal warming trends. By providing direct constraints on basin-scale averages with dense sampling in time, these data complement previous point measurements that alias local and transient temperature anomalies.


Introduction
The ocean is warming in response to accumulating greenhouse gases in the atmosphere, and its heat capacity dominates the climate system's thermal inertia.While the warming has been most pronounced in the surface ocean, heat transfer to the deep ocean importantly slows the climate change experienced at the surface (e.g., Hansen et al., 1981;Gregory, 2000;Held et al., 2010;Kostov et al., 2014).Roemmich et al. (2015), for example, estimated from Argo floats that the top 0.5 km of the global ocean warmed at a rate of 5 mK yr −1 between 2006 and 2013, whereas the layer between 0.5 and 2 km warmed at a rate of 2 mK yr −1 .Had the heat not been transferred below 0.5 km depth, the top 0.5 km would have warmed more than twice as rapidly (excluding feedbacks with the atmosphere and radiative transfer).It is thus important to constrain the heat transfer to the deep ocean.Because this transfer is achieved by processes that depend on uncertain parameterizations in climate models (the overturning circulation, mesoscale eddies, and diapycnal mixing), strong observational constraints are crucial.
It remains an observational challenge to isolate the small-amplitude and large-scale climate signal in the presence of much larger-amplitude but local fluctuations due to mesoscale eddies, internal waves, and other oceanic transients.To 2 km depth and since the mid-2000s, Argo floats have provided unprecedented coverage of the world ocean.Even with currently about 4000 floats, however, the Argo array still aliases mesoscale eddies, and regional estimates of warming rates remain uncertain (e.g., Wunsch, 2016;Dushaw, 2019;Fig. 1a).Additionally, and maybe more glaringly, Core Argo floats do not sample half of the ocean's volume: that below 2 km depth.While the Deep Argo array is expanding (e.g., Roemmich et al., 2019;Johnson et al., 2020), the extant record is short and limited to a few regions.Estimates of temperature change below 2 km rely heavily on sparse hydrographic sections that are sampled about once a decade (e.g., Roemmich and Wunsch, 1984;Purkey and Johnson, 2010;Desbruyères et al., 2016;Desbruyères et al., 2017;Volkov et al., 2017;Fig. 1a).State estimates constrain such changes by combining available observations with an ocean model (e.g., Wunsch and Heimbach, 2007), but the possibility of model error and the lack of an uncertainty estimate make it difficult to assess how reliable such estimates are.For example, Wunsch and Heimbach (2014) inferred widespread cooling of the abyssal ocean using the ECCO (v4r1) state estimate (relative to an assumed and corrected initial state), whereas direct estimates from repeat hydrographic sections tend to indicate a dominance of warming (e.g., Purkey and Johnson, 2010;Desbruyères et al., 2017).Clearly, better observational constraints are needed.
Recently, Wu et al. (2020) demonstrated that sound waves generated by earthquakes, so-called T waves, can be used to constrain basin-scale temperature change in the deep ocean.These measurements are based on the idea of "ocean acoustic tomography," which was originally proposed by Munk and Wunsch (1979) and successfully demonstrated at both planetary (Munk and Forbes, 1989) and basin scales (The ATOC Consortium, 1998).T waves propagate faster in warmer water, so changes in the average temperature encountered by these waves can be detected as changes in their travel time.The use of sound waves produced by earthquakes eliminates the need to deploy synthetic sound sources.Wu et al. used  the seismic station DGAR on Diego Garcia to receive T waves generated 2900 km away by earthquakes near Nias Island off Sumatra (Fig. 1a).To remove uncertainties in the source location and timing, repeating earthquakes were employed to extract the change in travel time between two events.A key advantage of such acoustic measurements over point measurements is that they intrinsically average the temperature change along the sound waves' path and therefore suffer much less from spatial aliasing.
In an accompanying manuscript (hereafter referred to as W23), we show that Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO) hydrophones are more sensitive T-wave receivers than land stations like DGAR, which allows the detection of smaller earthquakes at the CTBTO station H08 off Diego Garcia than is possible with DGAR data.This use of CTBTO data improves the time resolution of the inferred deep-ocean temperature change between Nias Island and Diego Garcia.In W23, we further use Nias Island T-waves received at the CTBTO station H01 off Cape Leeuwin to infer a time series of temperature change averaged along this 4400 km section extending into the extra-tropical ocean (Fig. 1a).Here, we make use of an additional advantage of the crisp arrivals of T waves in CTBTO hydrophone records: that changes in the travel time can be measured reliably at a number of different frequencies (Fig. 1b, Section 2).Because T waves at different frequencies are sensitive to different parts of the water column, this frequency dependence in the travel time change can be used to constrain the vertical structure of the temperature change (Fig. 2; cf., Munk and Wunsch, 1983;Shang, 1989).This approach is similar to that employed in remote sensing of the atmosphere (e.g., Fu et al., 2004), except that acoustic rather than electromagnetic waves are used.
A similar vertical-slice tomography scheme has been developed for acoustic measurements with synthetic sources (e.g., Munk andWunsch, 1979, 1982;Munk et al., 1995), but it has not been realized at a basin scale.The sources employed in the Acoustic Thermometry of Ocean Climate (ATOC) experiment had frequencies centered on 75 Hz and were effectively point sources at precisely known locations, allowing for a convenient eigenray description of the propagation from source to receiver.These rays have a known geometry and, under certain circumstances, can be resolved and identified in arrival patterns (Spiesberger et al., 1980;The Ocean Tomography Group, 1982;Cornuelle et al., 1993;The ATOC Consortium, 1998;Worcester et al., 1999).Relative changes in the arrival times of different rays can then be used to infer the vertical structure of temperature change, most easily in the range average.For example, the arrival time of steep rays that sample the surface mixed layer undergo a larger seasonal cycle than near-axial rays that are confined to the thermocline and deep ocean.On the bottom-mounted receivers used in ATOC, however, only relatively steep rays could be resolved and identified, and the information on the vertical structure of the temperature field was limited (Dushaw, 1999;Dushaw et al., 2009).
The propagation and arrival patterns of T waves are more complicated (e.g., Okal, 2008), but information on the vertical structure of the temperature change can still be extracted.Sensitivity kernels quantify how the arrival time changes in response to a temperature perturbation anywhere along the path and anywhere in the water column (Fig. 2; Wu et al., 2020).Consistent with expectations from simplified modal calculations (Fig. S1), these range-averaged sensitivity kernels shift upward in the water column for increasing frequencies (Fig. 2).A surface-intensified warming, for example, produces a larger reduction in travel time at high frequencies than at low frequencies.We use this frequency dependence to estimate a rough vertical structure of the range-averaged deep temperature anomalies.We find that the inferred structure of the anomalies along the two sections matches expectations based on the dynamics that produce the anomalies and is in general agreement with previous estimates, although the T-wave data tends to show stronger anomalies and more reliably captures perturbations due to mesoscale eddies (Section 3).

Inferring vertical structure
The starting point for the inference of vertical structure in the temperature anomalies are the time series of travel time anomalies at a few different low frequencies, constructed from a total of 11 479 repeating pairs that arise from 3457 earthquakes in 2005 to 2018 (Fig. 1b; W23).We choose 2.00 Hz, 3.00 Hz, and 4.00 Hz for the path to Diego Garcia and 2.50 Hz, 3.25 Hz, and 4.00 Hz for the path to Cape Leeuwin.A slightly higher minimum frequency is used for the latter path because 2.00 Hz T waves are less consistently received at Cape Leeuwin than at Diego Garcia.Measurements at higher frequencies are not reliable because the waveform correlation drops markedly beyond 4 Hz (e.g., Fig. 1b).For each frequency, we apply a Gaussian filter with width 0.5 Hz centered on that frequency before calculating the correlation function between the T-wave arrivals of an event pair, as described in W23.How these time series are obtained from measured travel time differences between repeating events is described in Text S1, and the cycle skipping corrections that are applied to the measurements are described in Text S2.
To turn the time series of travel time anomalies at different frequencies into an estimate of the evolving vertical structure of temperature anomalies, we perform a singular value decomposition (SVD) of the range-integrated sensitivity kernels (Fig. 2).The prob-lem is severely under-determined, so we can only hope to obtain a coarse estimate of the vertical temperature structure.Let K denote the matrix whose three rows contain the range-integrated sensitivity kernels at the three frequencies, discretized to a ∆z = 50 m grid.Then, at each event time t, we would like to invert τ(t) = KT(t)∆z for T(t), where τ(t) contains the T-wave travel time anomalies at the three frequencies, and T(t) contains the range-averaged temperature anomaly profiles on the same grid as the kernels.The SVD K = UΛV T yields the singular values λ 1 = 2.6 s K −1 km −1 , λ 2 = 0.21 s K −1 km −1 , and λ 3 = 0.022 s K −1 km −1 for the path to Diego Garcia and λ 1 = 4.6 s K −1 km −1 , λ 2 = 0.30 s K −1 km −1 , and λ 3 = 0.012 s K −1 km −1 for the path to Cape Leeuwin (see Text S1 for details).The rapid decay in the singular values is a result of the similarity of the sensitivity kernels at the chosen frequencies, i.e., their being nearly linearly dependent.Small singular values amplify errors (we only know the estimate τ, not the true τ), so a common trade-off must be made between resolution and precision.We choose to retain the first two singular vectors to obtain coarse vertical resolution with acceptable uncertainty: , where U 2 and V 2 consist of the first two columns of U and V, respectively, Λ 2 is the diagonal matrix containing λ 1 and λ 2 , and h = 5 km is a fixed reference depth.The estimate T2 (t) is then related to the true temperature field T(t) by T2 (t In the absence of errors in τ(t), the estimate T2 (t) is a projection of the true state T(t) onto the first two singular vectors, and the resolution matrix V 2 V 2 T determines to what degree features can be resolved by the available data.For both paths, only features between about 0.5 and 3 km depth can be captured, and the depth resolution is no better than about 1 km (Fig. 2).The projection coefficients are estimated from the data as c2 , with which the reconstructed temperature profile is the linear combination T2 (t) = V 2 c2 (t).(Note that the prior covariances for inferring τ(t) from the travel time changes between repeating events are chosen such that the three components of c(t) = h −1 Λ −1 U T τ(t) are uncorrelated.Any phase relations between them therefore arise entirely from the data.Details are given in Text S1.) Taking this SVD approach to inversion, we feign complete ignorance about the vertical structure of the range-averaged temperature anomalies.The inversion uses neither constraints from the physics that govern the temperature field nor prior information, say from previous measurements or model simulations.In particular, the inversion knows nothing of the typically strong surface intensification of temperature anomalies, which entails that inverted temperature anomalies tend to reach too deep and that upper-ocean anomalies can produce spurious oppositely signed deep anomalies (Fig. S6, S7).We nevertheless choose this agnostic approach to illustrate most simply what information is contained in the Twave data.Future work should combine this information from seismic data with prior information and other observations.
the path to Diego Garcia (Fig. 3a,b; 4a,b; Table 1).The first singular vectors have a similar shape as the kernels themselves, so their coefficients represent a weighted average of the deep temperature anomalies, with the weighting peaked at 1.5 km and 1.3 km for the paths to Diego Garcia and Cape Leeuwin, respectively (Fig. 2b,e).The second singular vectors have nodes at 1.7 km and 1.5 km for the two paths, so the projection onto them is more sensitive to where in the water column the anomalies are located.These projections inferred from T-wave travel time anomalies show general agreement with previous products.We compare our estimates to monthly interpolated Argo data (Roemmich and Gilson, 2009) and daily ECCO state estimate data (v4r4; Forget et al., 2015;ECCO Consortium, 2021).Using the sensitivity kernels, we infer travel time anomalies from the temperature anomalies of these products, interpolate those onto our event times, and subsequently treat them in the same way as the T-wave anomalies.
For the path to Diego Garcia, seasonal and sub-seasonal variations in both these projections tend to line up (Fig. 3a,b), indicating that the vertical structure inferred from the T-wave data is broadly consistent with these previous estimates.There are notable exceptions.The anomalies inferred from the T-wave data are stronger on average, particularly in the projection onto the first singular vector (Fig. 3a, Table 1).The projection onto the second singular vector (Fig. 3b) is more positive for the Argo product than for the T-wave data in the latter part of the time series, and it is less negative for ECCO than for the T-wave data in the early part of the time series, producing a stronger decadal trend in both previous products than inferred from the T-wave data (Table 1).
For the path to Cape Leeuwin, the seasonal variations in both projections are similar between those inferred from T-wave data and those from previous products (Fig. 4a,b).In contrast to the equatorial path to Diego Garcia, the annual signal here is much stronger than the semi-annual signal (Table 1).The T-wave data produces sizable spikes with a duration on the order of a month, which are typically missed by previous data (Fig. 4a,b; S7; S9).As discussed in W23, we interpret these spikes as resulting from mesoscale eddies traversing the path, most importantly those shed by the Leeuwin Current.The temperature anomalies of these eddies are largely confined to the thermocline (e.g., Fieux et al., 2005), which is consistent with the anomalies in the two projections appearing in phase.At times, these eddies happen to be captured by Argo floats (e.g., 2005-10 and 2010-10), but they are more typically missed, as expected from the Argo float coverage (Fig. 1).ECCO does not capture mesoscale eddies because it has too low a horizontal resolution.The overall stronger variability along this path to Cape Leeuwin implies that decadal trends are more uncertain (Table 1), although the uncertainties in the estimation from the three data sources are not independent.
Time series of the temperature reconstruction T2 illustrate the information on the vertical structure contained in the T-wave data.On the path to Diego Garcia, the inferred temperature anomalies are strongest in the upper 2 km but reach substantially below this depth (Fig. 3c).It cannot be inferred from the data whether these abyssal anomalies are real or whether they arise from the insufficient vertical resolution.(See Fig. S6, S7, S8, and S9 for projections of Argo and ECCO data onto the first two singular vectors.)Nevertheless, the anomalies display an upward phase propagation that is expected for long surface-generated equatorial waves that have a downward energy flux (e.g., Wunsch, 1977;Philander, 1978;McPhaden, 1982;Luyten and Roemmich, 1982;Reppin et al., 1999).This upward phase propagation is not apparent on the extratropical path to Cape Leeuwin (Fig. 4c), where the inferred anomalies are more strongly confined to above 2 km depth, consistent with mesoscale thermocline eddies.
In interpreting these results, it should be kept in mind that the sensitivity kernels are based on two-dimensional numerical simulations of the wave propagation from the source to the receiver (Wu et al., 2020).These simulations depend, albeit not sensitively, on assumptions about the thickness and properties of sediment layers and on the neglect of offgeodesic effects.There is therefore representational uncertainty in our estimates of temperature anomalies and their vertical structure arising from uncertainty in the kernels.More work is needed to better understand these effects.

Conclusions
The data presented here demonstrate that changes in T-wave travel times contain information on the vertical structure of the temperature anomalies encountered by these waves along their paths.This information can be extracted from travel time anomalies at a few low frequencies even though the sensitivity kernels at these frequencies have similar vertical structures (Fig. 2).
We here illustrated this vertical structure using a simple SVD inversion.In the future, the T-wave data should be combined with Argo and ship-based hydrographic data, either using a relatively simple mapping as typically employed for Argo data or using state estimation as in ECCO, which also allows one to incorporate additional constraints, for example from altimetry and gravitometry.T waves offer constraints on the large-scale temperature changes that are complementary to previous data.They intrinsically average in space, so they do not miss mesoscale eddies like the Argo array.They offer a dense sampling in time, which is important even for large-scale averages that still contain sizable anomalies induced by equatorial waves and mesoscale eddies.T waves offer constraints on the ocean below 2 km depth, which has been sparsely sampled in space and time by ship-based hydrographic surveys.The vertical resolution obtained from T-wave is relatively coarse, but they nevertheless constrain the vertical structure of the large-scale temperature anomalies.
In the present work, we restricted ourselves to low frequencies because the waveform correlation drops substantially at higher frequencies (Fig. 1b) and travel time changes cannot be extracted confidently.The T-wave signals have plenty of power at these higher frequencies, so noise is unlikely to be the problem.Instead, we speculate that the higher-frequency signals contain multiple vertical acoustic modes.The modes become more confined in the vertical as the frequency increases, so more modes escape interaction with the bottom and subsequent attenuation.Different modes experience different time shifts, so the waveform correlation drops if multiple modes contribute substantially.If this interpretation is correct, much more detailed information on the vertical structure could be extracted if we were able to separate the modes in the received signal, for example using a vertical hydrophone array (D'Spain et al., 2001).Deploying such arrays is routine, so future T-wave measurements could yield much stronger constraints on the vertical structure of the ocean's large-scale temperature variability.

Open research
The IMS hydrophone data are available directly from the CTBTO upon request and signing a confidentiality agreement to access the virtual Data Exploitation Centre (vDEC).All seismic data were downloaded through the IRIS Data Management Center (https://service.iris.edu/),including the seismic networks II (GSN; https://doi.org/10.7914/SN/II),MY, PS, and GE (https://doi.org/10.14470/TR560404).The Global Seismographic Network (GSN) is a cooperative scientific facility operated jointly by the Incorporated Research Institutions for Seismology (IRIS), the United States Geological Survey (USGS) and the National Science Foundation (NSF), under Cooperative Agreement EAR-1261681.The processing code is available at https://github.com/joernc/SOT.

Disclaimer
The views expressed in the paper are those of the authors and do not necessarily represent those of the CTBTO.

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 and origin-time corrections separately, which allows us to more naturally use multiple Twave 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.
Our primary goal is to infer time series of the travel time anomalies at the observed frequencies ω 1 , . . ., ω l , stacked into the vector These travel time anomalies, ascribed to changes in the ocean's temperature field, are differences between T-wave arrival anomalies and origin time corrections: where a (T)   i contains the T-wave travel time anomalies for each event time at the different 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 .
The travel time anomalies τ are related to the range-averaged temperature anomalies through the sensitivity kernel: 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 sensitivity matrix K = UΛV T .We normalize such that U T U = I 3 and h −1 V T V∆z = I l , where ∆z = 50 m is the vertical grid spacing and h = 5 km is a reference depth.With this normalization, 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 series 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: The pair matrices X (T) and X (L) take these differences between events at all T-wave and land stations involved.For example, if there were four pairs connecting five events observed at one T-wave section, we might have 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 The observed arrival time offsets, again measured relative to the cataloged event times, are then 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 : 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 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 appending their coefficients to the vector a and extending the matrices E, R, and D. For the design matrix, with where t contains the event times referenced to the midpoint between the first event and last event.For the prior covariance matrix, where Ξ is a diagonal matrix containing the prior variances for the trend or seasonal amplitudes 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, with B = X (T) t X (T) cos ωt X (T) sin ωt X (T) cos 2ωt X (T) sin 2ωt , 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 information required for a Gauss-Markov estimate: The desired estimate of the travel time anomalies is then τ = D ã, and the posterior covariance P can be propagated to this estimate as DPD T .Similarly, estimates of the projections onto singular vectors c, of the trends, and seasonal amplitudes of the signal can be obtained from ã, 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) 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 local currents, which contributes measurement error because the moorings are not navigated.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 mooring 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 of 30 m corresponds to a travel time anomaly of 0.02 s, which is equal to 2σ n .A more refined 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 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 propagating through Munk's canonical sound speed profile.The observed T-wave delays show a 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 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 with boundary conditions P = 0 at z = 0 (the surface) and dP/dz = 0 at z = −h (the bottom).We apply the normalization dz P 2 = h.The group slowness is 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 The change in group slowness, in contrast, is where (P 2 ) ω = 2PP ω , and P ω is defined by with the same boundary conditions as for P and the additional constraint that dz 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 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 temperature 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.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 sound speed profile, the sound speed perturbation is proportional to N 2 G.The dynamical mode G for vertical displacement is defined by with G = 0 at z = 0 and z = −h, and a normalization dz 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 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

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 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 surfaceintensified, 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 clusters with larger-magnitude low-frequency delays that suffer from forward or backward cycle skipping (Fig. S4, S5). 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 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 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).

Figure 1 :
Figure 1: Study region and example measurement of the frequency-dependent travel time change.(a) Locations of the Nias Island earthquakes and the hydrophone receivers at Diego Garcia and Cape Leeuwin.Also shown are the locations of all Argo profiles collected during the ten days following the 2005-03-28 Nias-Simeulue earthquake (orange dots) and the WOCE hydrographic sections (black dots).(b) Frequency-dependent waveform correlation function for the event pair 2005-04-12 09:56:39, 2008-10-11 12:32:32.The black dots indicate measurements for the central maximum as well as the two adjacent maxima that are used for cycle skipping correction (not needed for this pair).An origin time correction of 4.28 s was applied based on the waveforms received at land station PSI.

Figure 2 :
Figure 2: Sensitivity kernels and inference of vertical structure.Shown are (a), (d) the rangeaveraged sensitivity kernels at different frequencies, (b), (e) the first two singular vectors obtained from these range-averaged sensitivity kernels, and (c), (f) the resolution matrices of the SVD inversion employing the first two singular vectors.All these are shown for the path to Diego Garcia (top row) and that to Cape Leeuwin (bottom row).

Figure 3 :
Figure 3: Inferred temperature anomalies for the path to Diego Garcia.Shown are (a), (b) the projections onto the first two singular vectors c2 for the T-wave data and two reference products and (c) the SVD reconstruction of the range-averaged temperature anomaly profile T2 obtained from the T-wave data.Each dot in (a) and (b) corresponds to an estimate at an event time, and the shading indicates ±2σ confidence intervals around these estimates.Also shown are the estimated decadal trends of the projection coefficients c2 , including their ±2σ uncertainties.

Figure 4 :
Figure 4: Inferred temperature anomalies for the path to Cape Leeuwin.Shown are (a), (b) the projections onto the first two singular vectors c2 for the T-wave data and two reference products and (c) the SVD reconstruction of the range-averaged temperature anomaly profile T2 obtained from the T-wave data.Each dot in (a) and (b) corresponds to an estimate at an event time, and the shading indicates ±2σ confidence intervals around these estimates.Also shown are the estimated decadal trends of the projection coefficients c2 , including their ±2σ uncertainties.

Figure S2 :
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).

Figure S3 :Figure S6 :Figure S7 :
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 crosscorrelation 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 maximum (marked by triangle)

Figure S8 :
Figure S8: Temperature anomalies for the path to Diego Garcia projected onto the first two singular vectors.Shown are (a) the time series inferred from T-wave data, (b) the projections of anomalies from the Argo product, and (c) projections of anomalies from the ECCO product.

Figure S9 :
Figure S9: Temperature anomalies for the path to Cape Leeuwin projected onto the first two singular vectors.Shown are (a) the time series inferred from T-wave data, (b) the projections of anomalies from the Argo product, and (c) projections of anomalies from the ECCO product.
, due to noise affecting the cross-correlation