---
author:
- |
    Yanhua O. Yuan, Frederik J. Simons, and Jeroen Tromp\
    Department of Geosciences, Princeton University, Princeton, NJ 08544, USA\
    E-mail: yanhuay@princeton.edu
title: 'Double-difference adjoint seismic tomography'
---

We introduce a ‘double-difference’ method for the inversion for seismic wavespeed structure based on adjoint tomography. Differences between seismic observations and model predictions at individual stations may arise from factors other than structural heterogeneity, such as errors in the assumed source-time function, inaccurate timings, and systematic uncertainties. To alleviate the corresponding nonuniqueness in the inverse problem, we construct differential measurements between stations, thereby reducing the influence of the source signature and systematic errors. We minimize the discrepancy between observations and simulations in terms of the differential measurements made on station pairs. We show how to implement the double-difference concept in adjoint tomography, both theoretically and in practice. We compare the sensitivities of absolute and differential measurements. The former provide absolute information on structure along the ray paths between stations and sources, whereas the latter explain relative (and thus higher-resolution) structural variations in areas close to the stations. Whereas in conventional tomography a measurement made on a single earthquake-station pair provides very limited structural information, in double-difference tomography one earthquake can actually resolve significant details of the structure. The double-difference methodology can be incorporated into the usual adjoint tomography workflow by simply pairing up all conventional measurements; the computational cost of the necessary adjoint simulations is largely unaffected. Rather than adding to the computational burden, the inversion of double-difference measurements merely modifies the construction of the adjoint sources for data assimilation.\
**Key words:** Time-series analysis; Inverse theory; Tomography; Seismic tomography; Computational seismology; Wave propagation.

\[introduction\]I N T R O D U C T I O N
=======================================

The quality of tomographic inversions for seismic wavespeed structure is affected by uncertainties in our knowledge of the source-time function (the ‘source wavelet’ in exploration seismology), the (earthquake) source mechanism, location and origin time, inaccuracies in record timekeeping, network geometry, and systematic uncertainties [@Nolet2008]. All of these factors introduce intrinsic errors into synthetic modeling, even when accurate information on structure is available. Seismic tomography — in seeking to maximize the agreement between simulations and observations from all sources at all recording stations — is at risk of mapping errors of this kind into estimates of wavespeed variations [@Schaeffer+2013]. Jointly updating source terms and structural model parameters [@Pavlis+80; @Spencer+80; @Abers+91; @Thurber92], sequentially or iteratively, is a workable albeit expensive solution [@Widiyantoro+2000; @Panning+2006b; @Tian+2011; @Simmons+2012].

The fundamental ideas of making differential measurements were developed early on [e.g. @Brune+63; @Passier+95b]. The earthquake location community pioneered the concept of ‘double-difference’ inversions [@Poupinet+84; @Got+94], and a popular code, `hypoDD`, was introduced by [@Waldhauser+2000] to improve the accuracy of source locations. The double-difference tomography code `tomoDD` [@Zhang+2003b] produces high-resolution wavespeed models simultaneously with accurate event locations. Rather than considering pairs of earthquakes recorded at each station to reduce the effects of structural uncertainty on the source locations, differencing differential measurements between pairs of stations recording the same earthquake lessens the effect of uncertainties in the source terms on the determination of Earth structure [@Monteiller+2005; @Fang+2014]. Much as the double-difference technique has revolutionized earthquake (and other source) studies [e.g. @Rubin+99; @Rietbrock+2004; @Schaff+2004b], remaining the method of choice for high-resolution hypocenter determination, seismic tomography gains from fully embracing those concepts for structural inversions.

Until now, most of the double-difference approaches have operated within the confines of ray theory [but see @deVos+2013; @Maupin+2015]. In this paper, we bring double-difference inversions to full-waveform adjoint seismic tomography. Finite-frequency traveltime measurements drive a nonlinear inversion strategy that performs misfit gradient calculations via an elastic adjoint-state method. [@Dahlen+2000] and [@Hung+2000] approximated the Fréchet kernel for differential traveltimes, measured by cross-correlation of phases (e.g., *P*, *S*, *pP*, *SS*, ...) with ‘identical’ pulse shapes, by the difference of the individual single-phase kernels calculated asymptotically. [@Hung+2004] used the same approximation to measure the sensitivities of relative time delays between two nearby stations, which showed improved resolution of the upper-mantle structure beneath a regional array. [@Tromp+2005] brought the elastic adjoint method of [@Tarantola87b; @Tarantola88] into global seismology, calculating sensitivity kernels accurately and efficiently using a spectral-element method [@Komatitsch+2002a; @Komatitsch+2002b]. With these techniques it is now possible to incorporate essentially any type of seismic ‘measurement’ (absolute, relative, or differential) into global inversions for seismic wavespeed, boundary perturbations, and source parameters. In this paper, we derive a constructive theory for the incorporation of generic ‘double-difference’ measurements into adjoint seismic tomography.

The minimization of the difference of the difference between observations at distinct pairs of stations, and the difference between synthetics at the same pairs, over all station pairs, requires explicit mathematical expressions for the adjoint sources that are necessary for the numerical computation of the corresponding misfit gradient functionals. We derive those and show in numerical experiments that the gradients of the new misfit functional with respect to structural perturbations are relatively insensitive to an incorrect source signature and timing errors. These results stand in contrast to conventional tomography, which aims to minimize the difference between predictions and observations obtained from measurements made at individual stations.

A seismic phase observed at two ‘nearby’ (relative to the average source-receiver distance) stations will have ‘similar’ (using robust metrics) waveforms and ‘sensitivities’ (misfit gradients) to Earth structure. Station-relative differential measurements will reflect smaller-scale structural variations near to and in-between the station pairs. Structural inversions under the conventional formalism, unless actively ‘re’- or ‘pre’-conditioned to avoid doing so [@Curtis+97; @Spakman+2001; @Fichtner+2011b; @Luo+2015], tend to over-emphasize areas with high-density ray-path coverage relative to poorly sampled parts of the model. By the partial cancellation of common sensitivities, double-difference tomography, on the contrary, illuminates areas of the model domain where ray paths are not densely overlapping. Hence, using only one event (e.g., a teleseismic earthquake) recorded at a cluster of stations, conventional tomography can derive only very limited structural information, while double-difference tomography will resolve the structure in the instrumented area in greater detail, a benefit that accrues with the density of seismometer arrays [@Rost+2002; @Burdick+2014].

The computational cost for tomographic inversion is essentially unaffected by incorporating the double-difference concept into the conventional adjoint tomography workflow. All existing individual measurements are simply paired up, and only the construction of the adjoint sources requires modification. Data assimilation happens by back-propagating all the resulting adjoint sources simultaneously to compute the gradient kernels, exactly as for conventional adjoint tomography, without additional computational burden. Furthermore, the differential measurements themselves can be composed from any of the types used in conventional tomography. For example, they can be relative times calculated from a catalog of absolute arrival times [e.g., @VandeCar+90]; they may relative traveltime delays obtained from waveform cross-correlation analysis [e.g., @Luo+91], or from the cross-correlation of their envelopes [e.g., @Yuan+2015], and so on. The comparison need not be between the same phase observed at two distinct stations; it can involve different phases recorded on the same trace, or signals recorded at the same station but originating from two different events. In this paper, we use common waveform cross-correlation traveltimes as examples of incorporating differential measurements into adjoint-based tomography. For dispersive waves such as surface waves, differential frequency-dependent phase and amplitude measurements can be made using multitaper cross-spectral analysis. The details of that procedure are relegated to the Appendix.

Double-difference adjoint tomography relies on the assimilation of measurements, which we discuss how to group (e.g., via cluster analysis), and how to weight and stabilize (e.g., by regularization). We propose two approaches to cluster analysis and regularization. The first is based on source-station geometry, in which only station pairs whose spacing is comparable to the scale of the resolvable wavelengths are included in the double-difference data ensemble. The second is based on waveform similarity, as evaluated by the maximum normalized cross-correlation between two waveforms. When the waveform similarity exceeds a certain threshold, the pair of stations is integrated in the double-difference group, and the relative contributions to the misfit function from all qualifying station pairs are weighted by their similarity.

We demonstrate how to make double-difference measurements, and how to use them for adjoint tomography. Numerical experiments show the sensitivities of the double-difference data compared with conventional absolute measurements. We conduct tests with realistic network configurations, both on a global scale and at the scale of the North American continent, and for wavespeed structures that are either checkerboard synthetics or plausible Rayleigh-wave phase speed maps inspired by prior studies.

A D J O I N TT O M O G R A P H Y :T H EC L A S S I C A LA P P R O A C H
=======================================================================

We briefly review the principles of ‘conventional’ adjoint tomography using ‘absolute’ cross-correlation traveltimes. The material in this section is later used for comparison with the ‘double-difference’ adjoint tomography using ‘differential’ cross-correlation traveltimes.

The cross-correlation traveltime difference between a synthetic signal $s_i(t)$ and an observation $d_i(t)$ over a window of length $T$ is defined as $$\begin{aligned}
\label{eq:cc_sd}
  \Delta t_i&=\operatorname*{\arg\!\max}_\tau {\int_0^T}s_i(t+\tau) d_i(t) {\hspace{0.1em}\mathrm{d}t},
  \qquad
  \mbox{and we define}
  \qquad
  \gamma_{i}={\int_0^T}s_i(t+\tau) d_i(t) {\hspace{0.1em}\mathrm{d}t}.\end{aligned}$$ A positive $\Delta t_i$ indicates that the data waveform $d_i$ is advanced relative to the synthetic $s_i$, meaning that the wavespeed model that generated the synthetics is slower than the true model. The objective function to be minimized by ‘least-squares’ is the sum over all measurements of the squared traveltime shifts in eq. (\[eq:cc\_sd\]): $$\begin{aligned}
\label{eq:misfit_CC}
  \chi{_{\tiny \mathrm{cc}}}&={\frac{1}{2}}\sum_i  [\Delta t_i]^2
  .\end{aligned}$$ Explicit and complete derivations for the Fréchet derivatives of the terms in the cross-correlation traveltime misfit function in eq. (\[eq:misfit\_CC\]) were presented by various authors [e.g., @Luo+91; @Marquering+99; @Dahlen+2000]. Most central to the development, $$\label{dti}
  \delta\Delta t_i=
  \frac{{\int_0^T}{\partial}_t s_i(t)\,
    \delta s_i(t){\hspace{0.1em}\mathrm{d}t}}
       {{\int_0^T}{\partial}^2_t s_i(t)  s_i(t) {\hspace{0.1em}\mathrm{d}t}}
       .$$ To turn the expression for the traveltime perturbation, eq. (\[dti\]), into a useful expression for the perturbation of the misfit function in eq. (\[eq:misfit\_CC\]) requires a mechanism to relate a change in Earth properties (density and elastic constants or wavespeeds) to a change in the seismogram, $\delta s_i(t)$. Over the years, various formalisms were developed (acoustic, ray-based, modes-based, and numerical approaches), and they were amply discussed in the literature. In addition to the works cited above, we must still mention [@Zhao+2000], [@Chen+2007b] and [@Nissen-Meyer+2007a]. Here we follow [@Tromp+2005] in taking the numerical approach, which involves the action of an ‘adjoint source’. For an individual measurement of cross-correlation traveltime made at ${\mathbf{x}}_i$, the adjoint source $$\begin{aligned}
\label{eq:adj_CC}
  f_i^{\dagger}({\mathbf{x}},t)&= 
  \Delta t_i
  \frac{{\partial}_t s_i(T-t)}
       {{\int_0^T}{\partial}^2_t s_i(t) \, s_i(t) {\hspace{0.1em}\mathrm{d}t}} \, \delta ({\mathbf{x}}-{\mathbf{x}}_i)
       .\end{aligned}$$ Note the weighting by the traveltime anomalies $\Delta t_i$. The reverse-time synthetics generated at each of the stations from the corresponding adjoint sources in eq. ([\[eq:adj\_CC\]]{}) are summed and simultaneously back-propagated. The interaction of the adjoint with the forward-propagating wavefield then produces the gradients of the objective function in eq. (\[eq:misfit\_CC\]) with respect to the parameterized model perturbations $m$, leading to expressions that embody the essence of the inverse problem in seismic tomography [@Nolet96], namely $$\label{kp}
  \delta\chi{_{\tiny \mathrm{cc}}}=\int_\oplus K_m({\mathbf{x}}) \,m({\mathbf{x}}){\hspace{0.1em}\mathrm{d}^3{\mathbf{x}}},
  {\quad\mbox{where}\quad}K_m({\mathbf{x}})~\mbox{is the cross-correlation traveltime misfit sensitivity kernel for parameter}~m
  .$$ Positive cross-correlation traveltime differences $\Delta t_i$ from eq. (\[eq:cc\_sd\]) result in negative kernel values $K_m$ in eq. (\[kp\]). The model update required is in the direction opposite to $K_m$, which, for positive $\Delta t_i$ and negative $K_m$, implies that the model wavespeeds need to be sped up to reduce the misfit in eq. (\[eq:misfit\_CC\]). The integrations are carried out over the entire volume of the Earth, and the summation over the set of parameters $m$ is implied. The typical isotropic situation would be for eq. (\[kp\]) to involve the density and elastic moduli $\rho({\mathbf{x}})$, $\kappa({\mathbf{x}})$ and $\mu({\mathbf{x}})$, or the density and the elastic wavespeeds $\rho({\mathbf{x}})$, $\alpha({\mathbf{x}})$ and $\beta({\mathbf{x}})$, and so on. [@Tape+2007], [@Fichtner+2008], [@Zhu+2009], and others, describe in detail how the forward and adjoint wavefields interact to produce the various types of ‘misfit sensitivity kernels’ under different model parameterizations, and for the specific case considered in this section. An excellent recent overview is by [@Luo+2015].

A D J O I N TT O M O G R A P H Y :T H ED O U B L E - D I F F E R E N C EW A Y
=============================================================================

We use the term ‘double-difference’ measurement for the difference (between synthetics and observations) of differential measurements (either between pairs of stations or between pairs of events). From the new measurement we construct a new misfit function and derive its Fréchet derivatives and adjoint sources. In this section we focus on ‘measurements’ made by cross-correlation, and take ‘differential’ to mean ‘between pairs of stations, from a common source’. Our definitions can be relaxed later to apply to other types of measurements (e.g., waveform differences, envelope differences), and to refer more broadly to differential measurements, e.g., between different seismic phases or between different seismic sources observed at the same station.

Measurement
-----------

We consider the case of differential traveltimes calculated by cross-correlation of waveforms from a common source recorded at a pair of stations indexed $i$ and $j$. Let $s_i(t)$ and $s_j(t)$ denote a pair of synthetic waveforms, and let $d_i(t)$ and $d_j(t)$ be the corresponding pair of observations. The differential cross-correlation traveltimes, between stations, for the synthetic and the observation pairs are, respectively,

\[eq:cc\] $$\begin{aligned}
    \label{eq:cc_ss}
    {\Delta t{_{ij}}{^{\mathrm{syn}}}}&=\operatorname*{\arg\!\max}_\tau  {\int_0^T}s_i(t+\tau)s_j(t) {\hspace{0.1em}\mathrm{d}t},
    \qquad
    \mbox{and we define}
    \qquad
    \Gamma{_{ij}}(\tau)=  {\int_0^T}s_i(t+\tau)s_j(t) {\hspace{0.1em}\mathrm{d}t},
    \\
    \label{eq:cc_dd}
    {\Delta t{_{ij}}{^{\mathrm{obs}}}}&=\operatorname*{\arg\!\max}_\tau  {\int_0^T}d_i(t+\tau)d_j(t) {\hspace{0.1em}\mathrm{d}t},
    \qquad
    \mbox{and we define}
    \qquad
    \Lambda{_{ij}}(\tau)= {\int_0^T}d_i(t+\tau)d_j(t) {\hspace{0.1em}\mathrm{d}t}.
  \end{aligned}$$

A positive $\Delta t_{ij}$ indicates that the waveform recorded at station $j$ (i.e., $s_j$ or $d_j$) is advanced relative to the waveform recorded at station $i$ (i.e., $s_i$ or $d_i$). The difference of these differential traveltimes, between synthetics and observations, is the ‘double-difference’ traveltime measurement: $$\begin{aligned}
\label{eq:cc_DD}
  {\Delta\Delta}t{_{ij}}&= {\Delta t{_{ij}}{^{\mathrm{syn}}}}-{\Delta t{_{ij}}{^{\mathrm{obs}}}}.\end{aligned}$$ A positive ${\Delta\Delta}t{_{ij}}$ indicates that the advancement of $s_j$ over $s_i$ is larger than that of $d_j$ over $d_i$.

Misfit function
---------------

The misfit that we minimize is the sum of squares of the double-difference traveltime measurements, from eq. (\[eq:cc\_DD\]), between all station pairs: $$\begin{aligned}
\label{eq:misfit_CC_DD}
  \chi{_{\tiny \mathrm{cc}}}{^{\tiny \mathrm{dd}}}&={\frac{1}{2}}{\sum_{i}\sum_{j>i}}[{\Delta\Delta}t{_{ij}}]^2
  .\end{aligned}$$ The summation avoids double counting from the symmetry of the traveltime measurements, since $\Delta t{_{ji}}{^{\mathrm{syn}}}= -{\Delta t{_{ij}}{^{\mathrm{syn}}}}$ and $\Delta t{_{ji}}{^{\mathrm{obs}}}= -{\Delta t{_{ij}}{^{\mathrm{obs}}}}$.

Misfit gradient
---------------

The derivative of the differential objective function in eq. (\[eq:misfit\_CC\_DD\]) is $$\begin{aligned}
\label{eq:derivative_misfit_CC_DD}
  \delta \chi{_{\tiny \mathrm{cc}}}{^{\tiny \mathrm{dd}}}&= {\sum_{i}\sum_{j>i}}[{\Delta\Delta}t{_{ij}}]\,\delta{\Delta t{_{ij}}{^{\mathrm{syn}}}},\end{aligned}$$ where $\delta{\Delta t{_{ij}}{^{\mathrm{syn}}}}$ is the perturbation to the differential traveltime synthetic ${\Delta t{_{ij}}{^{\mathrm{syn}}}}$ caused by a model perturbation. As is customary we approximate the synthetic wavefield in the perturbed model to first order as the sum of the unperturbed wavefield $s(t)$ and a perturbed wavefield $\delta s(t)$, $$\label{eq:wavefield}
  \tilde{s}_i(t)=s_i(t)+\delta s_i(t) {\qquad\mbox{and}\qquad}\tilde{s}_j(t)=s_j(t)+\delta s_j(t).$$ Hence, from eqs (\[eq:cc\_ss\]) and (\[eq:wavefield\]), the cross-correlogram of the new seismograms $\tilde{s}_i(t)$ and $\tilde{s}_j(t)$ is, to first order in the perturbation, $$\begin{aligned}
\label{eq:Gamma}
  \tilde{\Gamma}{_{ij}}(\tau)
  &= {\int_0^T}\tilde{s}_i(t+\tau)\tilde{s}_j(t){\hspace{0.1em}\mathrm{d}t}\approx \Gamma{_{ij}}(\tau) + {\int_0^T}\delta s_i(t+\tau) s_j(t){\hspace{0.1em}\mathrm{d}t}+
  {\int_0^T}s_i(t+\tau) \delta s_j(t){\hspace{0.1em}\mathrm{d}t}= \Gamma{_{ij}}(\tau) + \delta \Gamma_i(\tau) + \delta\Gamma_j(\tau)
  .\end{aligned}$$ We introduced the notation $\delta \Gamma_i(\tau)$ for the cross-correlogram of $\delta s_i(t)$ an $s_j(t)$, and $\delta
\Gamma_j(\tau)$ for that between $s_i(t)$ and $\delta s_j(t)$, noting that $$\begin{aligned}
  \delta \Gamma_i(\tau) 
  &={\int_0^T}\delta s_i(t+\tau)s_j(t){\hspace{0.1em}\mathrm{d}t}= {\int_0^T}s_j(t-\tau) \delta s_i(t){\hspace{0.1em}\mathrm{d}t}{\qquad\mbox{and}\qquad}\delta  \Gamma_j(\tau) = {\int_0^T}s_i(t+\tau) \delta s_j(t){\hspace{0.1em}\mathrm{d}t}. 
  \label{eq:Gamma_part3}\end{aligned}$$ We recall from the definition of the differential cross-correlation traveltime in eq. (\[eq:cc\_ss\]) that the unperturbed cross-correlogram $\Gamma{_{ij}}(\tau)$ achieves its maximum at $\tau={\Delta t{_{ij}}{^{\mathrm{syn}}}}$, in other words, ${\partial}_\tau \Gamma{_{ij}}({\Delta t{_{ij}}{^{\mathrm{syn}}}}) = 0 $. If we expand the perturbed cross-correlogram $\tilde{\Gamma}{_{ij}}(\tau)$ in the vicinity of the unperturbed cross-correlation maximum differential time ${\Delta t{_{ij}}{^{\mathrm{syn}}}}$, keeping terms up to second order, we obtain $$\begin{aligned}
\label{eq:Gamma_expand}
  \tilde{\Gamma}{_{ij}}({\Delta t{_{ij}}{^{\mathrm{syn}}}}+\delta\tau)
  &= \Gamma{_{ij}}({\Delta t{_{ij}}{^{\mathrm{syn}}}})
  +\delta \Gamma_i({\Delta t{_{ij}}{^{\mathrm{syn}}}}) 
  +\delta \Gamma_j({\Delta t{_{ij}}{^{\mathrm{syn}}}}) 
  +\delta\tau \, {\partial}_\tau \delta \Gamma_i({\Delta t{_{ij}}{^{\mathrm{syn}}}}) 
  +\delta\tau \, {\partial}_\tau \delta \Gamma_j({\Delta t{_{ij}}{^{\mathrm{syn}}}})
  +{\frac{1}{2}}\delta\tau ^2 \: {\partial}^2_\tau \Gamma{_{ij}}({\Delta t{_{ij}}{^{\mathrm{syn}}}})
  .\end{aligned}$$ To find the perturbed time shift that maximizes the new cross-correlogram, we set its derivative with respect $\delta\tau$ to zero, thus requiring $$\begin{aligned}
\label{eq:Gamma_derivative}
  {\partial}_{\delta\tau}\tilde{\Gamma}{_{ij}}({\Delta t{_{ij}}{^{\mathrm{syn}}}}+\delta\tau) 
  &= \delta\tau \, {\partial}^2_\tau \Gamma{_{ij}}({\Delta t{_{ij}}{^{\mathrm{syn}}}})+
  {\partial}_\tau \delta \Gamma_i({\Delta t{_{ij}}{^{\mathrm{syn}}}}) +
  {\partial}_\tau \delta \Gamma_j({\Delta t{_{ij}}{^{\mathrm{syn}}}}) =0.\end{aligned}$$ The solution then yields the cross-correlation traveltime perturbation $\delta{\Delta t{_{ij}}{^{\mathrm{syn}}}}$ due to the model perturbation, namely $$\begin{aligned}
  \label{dtss}
  \delta{\Delta t{_{ij}}{^{\mathrm{syn}}}}&=-{\frac{\displaystyle{{\partial}_\tau \delta \Gamma_i({\Delta t{_{ij}}{^{\mathrm{syn}}}}) + {\partial}_\tau \delta \Gamma_j({\Delta t{_{ij}}{^{\mathrm{syn}}}})}}{\displaystyle{{\partial}^2_\tau \Gamma{_{ij}}({\Delta t{_{ij}}{^{\mathrm{syn}}}})}}} =  
      \frac{{\int_0^T}{\partial}_t s_j(t-{\Delta t{_{ij}}{^{\mathrm{syn}}}})\,\delta s_i(t){\hspace{0.1em}\mathrm{d}t}- {\int_0^T}{\partial}_t s_i(t+{\Delta t{_{ij}}{^{\mathrm{syn}}}})\,\delta s_j(t){\hspace{0.1em}\mathrm{d}t}} 
           {{\int_0^T}{\partial}^2_t s_i(t+{\Delta t{_{ij}}{^{\mathrm{syn}}}}) s_j(t){\hspace{0.1em}\mathrm{d}t}}.\end{aligned}$$ We introduce a notation for the denominator, $$N{_{ij}}={\int_0^T}{\partial}^2_t s_i(t+{\Delta t{_{ij}}{^{\mathrm{syn}}}}) s_j(t){\hspace{0.1em}\mathrm{d}t},$$ and rewrite eq. (\[eq:derivative\_misfit\_CC\_DD\]), the derivative of the differential cross-correlation objective function in eq. (\[eq:misfit\_CC\_DD\]), as $$\begin{aligned}
\label{dchi}
  \delta \chi{_{\tiny \mathrm{cc}}}{^{\tiny \mathrm{dd}}}&= {\int_0^T}\left\{
  \sum_{i}  \left[ \sum_{j>i}   
    {\frac{{\Delta\Delta}t{_{ij}}}{N{_{ij}}}}{\partial}_t s_j(t-{\Delta t{_{ij}}{^{\mathrm{syn}}}}) \right] \delta s_i(t)   
  -\sum_{j}  \left[  \sum_{i<j}
    {\frac{{\Delta\Delta}t{_{ij}}}{N{_{ij}}}}{\partial}_t s_i(t+{\Delta t{_{ij}}{^{\mathrm{syn}}}})
    \right]  \delta s_j(t)  \right\}{\mathrm{d}t}.\end{aligned}$$ Any scheme by which we minimize eq. (\[eq:misfit\_CC\_DD\]) proceeds until model perturbations no longer produce meaningful adjustments $\delta{\Delta t{_{ij}}{^{\mathrm{syn}}}}$ in the differential measurement on the synthetics ${\Delta t{_{ij}}{^{\mathrm{syn}}}}$, in which case $\delta \chi{_{\tiny \mathrm{cc}}}{^{\tiny \mathrm{dd}}}$ effectively vanishes and a minimum $\chi{_{\tiny \mathrm{cc}}}{^{\tiny \mathrm{dd}}}$ is reached.

Adjoint source
--------------

The translation of the gradient expression (\[dchi\]) into an algorithm for which the spectral-element forward-modeling software can be used to produce linearized expressions of the kind in eq. (\[kp\]) again draws upon the Born approximation and the reciprocity of the Green functions of the wave equation [@Tromp+2005]. An alternative derivation is through the lens of wave-equation constrained functional minimization [@Liu+2006]. Either route leads to a paired set of adjoint sources, formulated with respect to the synthetic wavefields $s_i$ and $s_j$ recorded at locations ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$. Considering all possible such pairs $i$ and $j$, we write the adjoint sources explicitly as

\[eq:adj\_CC\_DD\] $$\begin{aligned}
    \label{eq:adj_CC_DD_1}  
    f_i^{\dagger}({\mathbf{x}},t)
    &=+\sum_{j>i} {\frac{{\Delta\Delta}t{_{ij}}}{N{_{ij}}}}\, {\partial}_t s_j \big( T-[t-{\Delta t{_{ij}}{^{\mathrm{syn}}}}] \big)  
    \delta ({\mathbf{x}}-{\mathbf{x}}_i)
    ,\\
    \label{eq:adj_CC_DD_2}
    f_j^{\dagger}({\mathbf{x}},t)
    &=-\sum_{i<j} {\frac{{\Delta\Delta}t{_{ij}}}{N{_{ij}}}}\, {\partial}_t s_i \big( T- [t+{\Delta t{_{ij}}{^{\mathrm{syn}}}}] \big) 
    \delta ({\mathbf{x}}-{\mathbf{x}}_j)
    .
  \end{aligned}$$

The comparison of the expression for the adjoint sources of the classical case in eq. (\[eq:adj\_CC\]) with its double-difference counterpart in eqs (\[eq:adj\_CC\_DD\]) reveals that the former only involves one waveform per station, whereas in the latter case the adjoint source for each station comprises the waveforms of all the other stations with which it is being compared. As in the classical case, however, the final adjoint source for double-difference measurements is the sum of the contributions from all stations, which are thus back-propagated simultaneously to obtain the adjoint kernels of the objective function in eqs (\[eq:misfit\_CC\_DD\]). All the pairs of differential measurements are included in the adjoint process, which is otherwise identical to that used in conventional adjoint tomography. Since we linearly combine adjoint sources rather than combining individual adjoint kernels for each data ensemble, only one adjoint simulation is required to numerically evaluate the gradient of the double-difference misfit function, and once again we obtain an expression of the type in eq. (\[kp\]), with the summation over the parameter type $m$ implied: $$\label{kp2}
  \delta\chi{_{\tiny \mathrm{cc}}}{^{\tiny \mathrm{dd}}}=\int_\oplus K{^{\tiny \mathrm{dd}}}_m({\mathbf{x}}) \,m({\mathbf{x}}){\hspace{0.1em}\mathrm{d}^3{\mathbf{x}}},
  {\quad\mbox{with}\quad}K{^{\tiny \mathrm{dd}}}_m({\mathbf{x}})~\mbox{the double-difference cross-correlation traveltime misfit sensitivity kernel for}~m
  .$$ A positive ${\Delta\Delta}t{_{ij}}$ defined in eq. (\[eq:cc\_DD\]) results in positive kernel values (\[kp2\]) along the path between the source and station $j$, and negative kernel values along the path towards station $i$. Therefore, to decrease the double-difference misfit function in eq. (\[eq:misfit\_CC\_DD\]), the model update moves in the negative kernel direction, increasing the wavespeed to station $i$, reducing it towards station $j$, and decreasing the advancement of $s_j$ over $s_i$.

Cluster analysis and regularization
-----------------------------------

From $N$ individual measurements a total of ${{N}\choose{2}}=n(n-1)/2$ double-difference measurements can be created. To enhance the stability of the inverse problem and reduce the computational cost of the data pairings, we can regularize the problem via cluster analysis in which only pairs that are well linked to each other are included. In general this can be achieved by incorporating pairs of stations that are relatively ‘close’, which enhances accuracy and improves tomographic resolution. The maximum separation, however, depends on the wavelength $\lambda$ and the path length $L$ of the phases considered, the scale length of the heterogeneity, and the overall geometry of the illumination of the two stations by seismic sources. In practice a physical distance proportional to $\sqrt{\lambda L}$, the width of the first Fresnel zone, is recommended [e.g., @Woodward92; @Baig+2003]. We can furthermore exclude any station pairs whose spatial offset is much smaller than the expected scale lengths of the structural variations.

Another possible method of clustering is to select station pairs with ‘similar’ signals, as measured by the correlation $$\begin{aligned}
\label{eq:cc_coherency}
  r{_{ij}}&=
  \frac{\Lambda{_{ij}}({\Delta t{_{ij}}{^{\mathrm{obs}}}})}{\sqrt{\Lambda_{ii}(0) \Lambda_{jj}(0)}},
  \qquad
  -1\le r{_{ij}}\le 1
  ,\end{aligned}$$ where we recall the definition of $\Lambda_{ij}$ from eq. (\[eq:cc\_dd\]). Similarity pairs can be defined on the basis of threshold values for the correlation [@VandeCar+90], or the squared correlation can be used as weights [@Waldhauser+2000]. The similarity criterion can be used alone, or in combination with the relative physical proximity of the stations.

T H E A D V A N T A G EO FD I F F E R E N T I A LM E A S U R E M E N T S
========================================================================

In this section we briefly list some of the properties of differential traveltime measurements that make them good candidates for being less affected by source uncertainties than traditional absolute measurements. In the next section, these properties will be verified and illustrated by numerical experiments.

Scaling invariance
------------------

The cross-correlation measurement in eq. (\[eq:cc\_ss\]) is invariant to scaling of the seismograms $s_i(t)$ and $s_j(t)$ by a constant factor $\alpha$: $$\begin{aligned}
\label{eq:cc_scaling}
  \operatorname*{\arg\!\max}_\tau {\int_0^T}[\alpha s_i(t+\tau)][\alpha s_j(t)] {\hspace{0.1em}\mathrm{d}t}=\operatorname*{\arg\!\max}_\tau  \left[\alpha^2 {\int_0^T}s_i(t+\tau)  s_j(t) {\hspace{0.1em}\mathrm{d}t}\right]
  = \operatorname*{\arg\!\max}_\tau {\int_0^T}s_i(t+\tau) s_j(t) {\hspace{0.1em}\mathrm{d}t}={\Delta t{_{ij}}{^{\mathrm{syn}}}}.\end{aligned}$$ The scale invariance directly removes any sensitivity to incorrect seismic moment estimates.

Shift invariance
----------------

The time shift calculated by maximizing the cross-correlogram in eq. (\[eq:cc\_ss\]) is unchanged if time is globally shifted by some constant $\tau_0$: $$\begin{aligned}
\label{eq:cc_shift}
  \operatorname*{\arg\!\max}_\tau {\int_0^T}s_i(t+\tau-\tau_0) \: s_j(t-\tau_0) {\hspace{0.1em}\mathrm{d}t}= \operatorname*{\arg\!\max}_\tau {\int_0^T}s_i(t+\tau)\; s_j(t) {\hspace{0.1em}\mathrm{d}t}={\Delta t{_{ij}}{^{\mathrm{syn}}}}.\end{aligned}$$ The shift invariance directly removes any sensitivity to errors in origin time.

Source wavelet invariance
-------------------------

Knowledge of the source-time-function or ‘wavelet’ is of course crucial to making accurate measurements. Differential measurements, from the same seismic event to a pair of stations, are more robust to errors in the source wavelet than absolute ones. They should also be relatively robust to errors in the source’s focal mechanism, unless incorrect assumptions should place one of the event-station paths — but not both — on the wrong side of the nodal plane. The cross-correlation time delays between two synthetic seismograms calculated with the wrong source information should be relatively robust as long as we can view the ‘wrong’ records as filtered, or otherwise mostly linearly mapped, from what would be the ‘right’ seismograms. As long as the phases targeted by the cross-correlation time window remain well separated temporally from neighboring phases, in either case, they could be remapped to the same impulsive arrivals after deconvolution of their respective input source wavelets, hence preserving the differential delay time.

Differential sensitivity versus sensitivity of differential measurement
-----------------------------------------------------------------------

Eq. (\[dtss\]) contains two terms that contribute to the traveltime sensitivity, the Fréchet derivatives for the differential cross-correlation traveltime. Both of these mix contributions from both seismograms $s_i$ and $s_j$. Now suppose waveforms of $s_i(t)$ and $s_j(t)$ are identical upon time shifting, i.e., they are of the ‘identical pulse shape’ type (perhaps achieved by pre-processing, e.g., to equate the Maslov indices of *P* and *PP*, or *S* and *SS*) as identified by [@Dahlen+2000]. In that case $$\begin{aligned}
  s_i(t) \approx s_j(t-{\Delta t{_{ij}}{^{\mathrm{syn}}}}) {\qquad\mbox{and}\qquad}s_j(t) \approx s_i(t+{\Delta t{_{ij}}{^{\mathrm{syn}}}})
  ,\end{aligned}$$ and in that case eq. (\[dtss\]) would reduce to a form that does not mix station indices, and we learn by comparison with eq. (\[dti\]) that $$\begin{aligned}
  \label{eq:derivative_s1_app}
  \delta{\Delta t{_{ij}}{^{\mathrm{syn}}}}=  
  \frac{{\int_0^T}{\partial}_t s_i(t)\delta s_i(t){\hspace{0.1em}\mathrm{d}t}}{{\int_0^T}{\partial}^2_t s_i(t)s_i(t){\hspace{0.1em}\mathrm{d}t}}
  -
  \frac{{\int_0^T}{\partial}_t s_j(t)\delta s_j(t){\hspace{0.1em}\mathrm{d}t}}{{\int_0^T}{\partial}^2_t
    s_j(t)s_j(t){\hspace{0.1em}\mathrm{d}t}}
  =
  \delta\Delta t_i
  -\delta\Delta t_j
  . \end{aligned}$$ In other words, in the case of ‘identical’ waveforms, the Fréchet sensitivity kernel for their differential traveltime can be approximated by taking the difference of the individual Fréchet sensitivity kernels for each of the absolute traveltimes, and the classical formalism of adjoint tomography can be used to perform double-difference adjoint tomography, albeit approximately. This is the approach adopted by, among others, [@Hung+2004].

Conditional statements of the kind that validated the applicability of eq. (\[eq:derivative\_s1\_app\]) for using differential phases (either for different phases at the same station, or from the same phase at different stations) were made by [@Dahlen+2000] (compare our eq. \[dtss\] with their eq. 75). These assumptions are, however, not generally valid. In contrast, our eqs (\[eq:adj\_CC\_DD\]) are proper generalizations of the differential traveltime adjoint sources of [@Tromp+2005] (their eq. 61) which do not suffer similar such restrictions — save for the general adequacy of the Born approximation, which appears practically unavoidable though perhaps not theoretically unassailable [@Hudson+81; @Panning+2009; @Wu+2014b].

As a further point of note, it is computationally more efficient and algorithmically less complicated to evaluate the double-difference adjoint kernel using eqs (\[eq:adj\_CC\_DD\]) than via eqs (\[eq:adj\_CC\]) and (\[eq:derivative\_s1\_app\]). For example, for a total number of $N$ measurements, one only needs one adjoint simulation to numerically calculate the gradient of the differential misfit under ${{N}\choose{2}}$ possible combinations. In contrast, by using the approximation in eq. (\[eq:derivative\_s1\_app\]), one needs to compute the Fréchet sensitivity kernel for each individual measurement, which requires $N$ adjoint simulations in total. Subsequently these sensitivity kernels need to be combined by weighted with the corresponding double-difference traveltime measurements, and summed up for all grouped kernels. In conclusion, the misfit sensitivity for differential measurements for use in double-difference adjoint tomography is best obtained, both theoretically and in practice, via the paired adjoint sources of eqs (\[eq:adj\_CC\_DD\]).

N U M E R I C A LE X P E R I M E N T S
======================================

Three types of validation experiments are conducted here. We consider input checkerboard patterns but also a realistic phase speed distribution, for sensor configurations that proceed in complexity from two stations to simple arrays, and then to realistic array configurations.

Checkerboard model tests, local networks
----------------------------------------

We first conduct a numerical experiment with strong lateral heterogeneity, with a 2-D ‘checkerboard’ input model as shown in Fig. \[fig:Tape\_test1\_model\]. We use membrane surface waves [@Tanimoto90; @Peter+2007] as an analogue for short-period (10–20 s) surface waves. We use `SPECFEM2D` [@Komatitsch+98] to model the wavefield. Absorbing boundary conditions are employed on all sides.

Experiment I, as shown in Fig. \[fig:Tape\_test1\_model\], involves one source (filled star), and two stations (open circles labeled 1 and 2) that are aligned at different epicentral distances along a common source-receiver ray path. The target model of interest is a 480 km $\times$ 480 km checkerboard expressed as four circular *S*-wave speed anomalies of alternating sign, shown in Fig. \[fig:Tape\_test1\_model\]a. The initial model is a featureless homogeneous model with *S*-wave speed of 3500 m/s, shown in Fig. \[fig:Tape\_test1\_model\]b. We uniformly discretize each dimension using 40 elements of individual length 12 km. The source wavelet used is the first derivative of a Gaussian function with a dominant period of 12 s. The sampling rate is 0.06 s and the total recording length is 4.8 s. The exact source-time function and event origin time were used in the computation.

Fig. \[fig:Tape\_test1\_data\]a and Fig. \[fig:Tape\_test1\_data\]b show the displacement seismograms recorded at stations 1 and 2, respectively. The black traces marked ‘data’ are calculated in the checkerboard model of Fig. \[fig:Tape\_test1\_model\]a. The green traces labeled ‘synthetic’ are predictions in the homogeneous model of Fig. \[fig:Tape\_test1\_model\]b. We evaluate the time shift between the predictions (green) and the observations (black) at each of the stations via the conventional cross-correlation approach, as in eq. (\[eq:cc\_sd\]). These measurements are $\Delta t_1=4.08$ s ($d_1$ is advanced 4.08 s relative to $s_1$) and $\Delta t_2=6.60$ s ($d_2$ is advanced 6.60 s relative to $s_2$). Following eq. (\[eq:misfit\_CC\]) in this example, the objective function is the sum of the squares of the individual time shifts $\chi{_{\tiny \mathrm{cc}}}={\frac{1}{2}}([\Delta t_1]^2+[\Delta
  t_2]^2)=30.1$ s$^2$. Fig. \[fig:Tape\_test1\_data\]c and Fig. \[fig:Tape\_test1\_data\]d show the adjoint sources calculated from eq. (\[eq:adj\_CC\]), in non time-reversed coordinates, as blue lines labeled ‘cc’.

The double-difference cross-correlation approach measures, via eq (\[eq:cc\_ss\]), the time shift between the predictions made at station 1 and 2 (green traces in Figs \[fig:Tape\_test1\_data\]a–b) and, via eq. (\[eq:cc\_dd\]), between the observations at those same stations (black lines in Figs \[fig:Tape\_test1\_data\]a–b). The first of these is ${\Delta t_{12}{^{\mathrm{syn}}}}=-32.34$ s ($s_2$ is delayed 32.34 s relative to $s_1$), and the second ${\Delta t_{12}{^{\mathrm{obs}}}}=-29.76$ s ($d_2$ is delayed 29.76 s relative to $d_1$). Their ‘double’ difference, from eq. (\[eq:cc\_DD\]), yields ${\Delta\Delta}t_{12}=-2.58$ s (the delay of $s_2$ over $s_1$ is 2.58 s larger than that of $d_2$ over $d_1$). The double-difference objective function of eq. (\[eq:cc\_DD\]) is the sum of the squared double differences over all available pairs, which in this case simply evaluates to $\chi{_{\tiny \mathrm{cc}}}{^{\tiny \mathrm{dd}}}={\frac{1}{2}}[{\Delta\Delta}t_{12}]^2=
3.33$ s$^2$. Fig. \[fig:Tape\_test1\_data\]c and Fig. \[fig:Tape\_test1\_data\]d show the adjoint sources calculated from eqs (\[eq:adj\_CC\_DD\_1\]) and (\[eq:adj\_CC\_DD\_2\]), again in natural time coordinates, as red lines, which are labeled ‘cc\_dd’.

Fig. \[fig:Tape\_test1\_kernel\] shows the misfit sensitivity kernels for both types of measurements. The cross-correlation traveltime *S*-wave speed adjoint misfit kernel $K_\beta({\mathbf{x}})$ shown in Fig. \[fig:Tape\_test1\_kernel\]a is most sensitive to the area of effective overlap of the ray paths from the source to either station. This behavior is indeed explained by the nature of the kernel, which is obtained by summation of the individual adjoint kernels for each station. The desired model update is in the opposite kernel direction, increasing the wavespeed to station 1 and to station 2. In contrast, the double-difference cross-correlation traveltime *S*-wave speed adjoint misfit kernel $K{^{\tiny \mathrm{dd}}}_\beta({\mathbf{x}})$ shown in Fig. \[fig:Tape\_test1\_kernel\]b displays a dominant sensitivity in the area in-between the two stations. The ‘double’ differencing of the traveltime metrics between the distinct stations, cancels out much of the sensitivity common to both ray paths, the area between the source and station 1 in this case. The traveltime delay of the waves recorded at station 2, relative to station 1, mainly reflects the average *S*-wave speed between stations 1 and 2. The finite-frequency character of the measurements remains apparent from the ‘swallowtails’ visible in the area between the source and station 1. The opposite kernel direction required for model update is to slow down along the path from the source to station 1, and to speed up towards station 2. Since the path of station 1 overlaps that of station 2, the resulting effect is to increase the wavespeed between station 1 and station 2.

Experiment II, as shown in Fig. \[fig:Tape\_group2\_model\], again involves one source (filled star) and two stations (open circles), but now the latter have been arrange to lie at the same epicentral distance but at different azimuthal directions relative to the source.

In Experiment IIA we consider the ideal situation when source information is known. The exact timing and the exact source wavelet (again, the first derivative of a Gaussian with dominant period of 12 s) are used for synthetic modeling. Fig. \[fig:Tape\_test4\_data\]a and Fig. \[fig:Tape\_test4\_data\]b show the displacement seismograms recorded at both stations, with the checkerboard-model traces in black and the homogeneous-model traces in green. The cross-correlation time shift measurements are now $\Delta t_1=0.54$ s ($d_1$ is advanced 0.54 s relative to $s_1$) and and $\Delta t_2=-0.54$ s ($d_2$ is delayed 0.54 s relative to $s_2$). These are small and almost identical as the ray paths sample almost equal path lengths of slow and fast wavespeed anomalies. Fig. \[fig:Tape\_test4\_data\]c and Fig. \[fig:Tape\_test4\_data\]d show the adjoint sources for this model setup. In blue, for the traditional cross-correlation time shift between predictions and observations, in red for the difference of the differential cross-correlation time shifts between stations, where ${\Delta t_{12}{^{\mathrm{syn}}}}=0.00$ s (as $s_1$ and $s_2$ arrive at the same time), ${\Delta t_{12}{^{\mathrm{obs}}}}=-1.08$ s (as $d_2$ is delayed by 1.08 s relative to $d_1$), and the ‘double’ difference ${\Delta\Delta}t_{12}=1.08$ s. Fig. \[fig:Tape\_test4\_kernel\] shows the misfit sensitivity kernels for both types of adjoint measurements.

![image](Fig/Tape_test1_case1_adj){width="\traces"}

![image](Fig/Tape_test2_case1_adj){width="\traces"}

![image](Fig/Tape_test2_case2_adj){width="\traces"}

The $K_m({\mathbf{x}})$ and $K_m{^{\tiny \mathrm{dd}}}({\mathbf{x}})$ kernels, at first glance, look very similar, but their physical meanings are very different. The traditional measurements provide absolute information on wavespeed structure along the ray paths, whereas the double-difference measurements are sensitive to relative structural variations. From the classical prediction-observation cross-correlation measurements, we learn that the surface waves predicted to be recorded at station 1 arrive 0.54 s later than the observations, and those predicted at station 2 arrive 0.54 s earlier than the observations. To reduce the misfit function (in eq. \[eq:misfit\_CC\]), the wavespeed along the ray path from the source to station 1 should be increased, and the wavespeed between the source and station 2 decreased. On the other hand, from the differential cross-station traveltime measurements, we deduce that there is a 1.08 s advance in the arrival of the observed waves at station 1 relative to station 2, but the synthetic records arrive at exactly the same time at both stations, based on the lag time of their cross-correlation maximum. Of course here too, to reduce the misfit (in eq. \[eq:misfit\_CC\_DD\]), the wavespeed from the common source to station 1 should be made faster than that to station 2.

In Experiment IIB we treat the more realistic scenario where an incorrect source wavelet is used for the synthetic modeling. Fig. \[fig:Tape\_test5\_data\]a and \[fig:Tape\_test5\_data\]b show the displacement synthetics. The observations (black) are in the checkerboard model and use the first derivative of the Gaussian function for a source wavelet. The predictions (green) use the second derivative of the Gaussian, the ‘Ricker’ wavelet. Both source wavelets share the same excitation time and dominant period of 12 s. The cross-correlation time shifts between the synthetics and observations are $\Delta t_1=-2.64$ s and $\Delta
t_2=-3.72$ s. These are markedly different from the measurements obtained above using the correct source wavelet (recall those were $\Delta t_1=0.54$ s, $\Delta t_2=-0.54$ s). In contrast, the double-difference time shift between station 1 and station 2 is unchanged from that which uses the correct source wavelet, ${\Delta\Delta}t_{12}= 1.08$ s. Fig. \[fig:Tape\_test5\_kernel\] shows the misfit kernels.

The effects of an incorrect source signature on the traditional cross-correlation measurements can be clearly seen by comparing the adjoint kernel in Fig. \[fig:Tape\_test5\_kernel\]a with that shown in Fig. \[fig:Tape\_test4\_kernel\]a. These are markedly different. Fortunately, the double-difference measurement approach yields relatively consistent adjoint kernels, based on the comparison of the kernel shown in in Fig. \[fig:Tape\_test5\_kernel\]b with the one using exact source signature in Fig. \[fig:Tape\_test4\_kernel\]b. The latter two are nearly identical, which is indicative of desirable robustness in the inversion.

In Experiment IIC we next consider the scenario of having an incorrect origin time, an advance of 0.96 s expressed by the synthetics. Fig. \[fig:Tape\_test6\_data\]a and Fig. \[fig:Tape\_test6\_data\]b show the records, with the same color conventions as before. It comes as no surprise that the cross-correlation time shifts between predictions and observations are biased. We obtain $\Delta t_1=-0.42$ s and $\Delta t_2=-1.50$ s, a spurious shift of 0.96 s for each measurement. The differential cross-correlation time shifts between stations, however, are not affected at all, since the predictions at both stations, wrong as they both are, are shifted in exactly the same way with respect to the true source origin time. Fig. \[fig:Tape\_test6\_kernel\] shows the misfit kernels. The conventional cross-correlation kernel, in Fig. \[fig:Tape\_test6\_kernel\]a, are clearly affected by the timing error, while the double-difference ones, in Fig. \[fig:Tape\_test6\_kernel\]b, are not biased at all by any such inaccuracies.

In Experiment III, we consider one earthquake and an array of stations inspired by [@Tape+2007], using location information from 132 broadband receivers in the Southern California Seismic Network (SCSN). The target and initial models and the station-receiver geometry are shown in Fig. \[fig:Tape\_group3\_model\]. An iterative inversion based on conventional cross-correlation misfit measurements and eq. (\[kp\]) results in the estimated structure shown in in Fig. \[fig:Tape\_test7\_model\]a. An inversion based on double-difference cross-correlation misfit measurements and eq. (\[kp2\]) is shown in Fig. \[fig:Tape\_test7\_model\]b. The final models show great differences between them. The misfit gradient obtained with the conventional approach provides very limited structural information, and the ‘final’ model displays prominent ‘streaks’ between the stations and the source. In contrast, the double-difference method is largely free from such artifacts and yields significant detail on structural variations in the area under consideration.

Our checkerboard synthetic experiments were designed to demonstrate the ease and power of the double-difference concept in adjoint tomography, from making the measurements to calculating the kernels, to inverting for the best-fitting models. It is immediately apparent that the double-difference approach provides powerful interstation constraints on seismic structure, and is less sensitive to error or uncertainty in the source. Even using just a single earthquake, the double-difference method enables the iterative spectral-element adjoint-based inversion procedure to capture essential details of the structure to be imaged.

In experiments (not shown here) where the starting model is homogeneous but does not have the same mean wavespeed as the target model, the double-difference approach manages to recover much of the absolute and relative structure interior to the array, although regions outside of the region encompassing the bulk of the stations remain unupdated. In such cases, a hybrid method that combines both traditional and double-difference measurements brings improvements throughout the imaged region.

![image](Fig/Tape_test2_case3_adj){width="\traces"}

Realistic phase speed input models, global networks
---------------------------------------------------

To test realistic length scales in a global tomography model, we use the global phase velocity map of [@Trampert+95] at periods between 40 s and 150 s, which we mapped to a simplified, unrealistic geometry of a plane with absorbing boundaries. In Experiment IV we model membrane surface waves in this model, as shown in Fig. \[global\_test1\]a. The 32 selected earthquakes are shown by red stars, and white circles depict 293 station locations from the global networks.

We use the `SPECFEM2D` code to model the wavefield. We use 40 elements in latitude and 80 elements in longitude with an average of 4.5$^\circ$ in element size in each direction. The rescaling corresponds to dimensions of one meter per latitudinal degree ($^\circ$) of the original. The source wavelet is a Ricker function with 260 Hz dominant frequency, which should be able to resolve structure down to a scale length of approximately 6.6$^\circ$. The sampling rate is 26 $\mu$s and the recording length 125 ms.

The initial model is homogeneous with a phase speed of 4500 m/s. Using single-station cross-correlation traveltime measurements, the inversion yields a model that has very limited resolution, shown in Fig. \[global\_test1\]b. From double-difference cross-correlation measurements between stations, the inversion reveals structural variations with much greater detail, as shown in Fig. \[global\_test1\]c. The conventional tomographic method appears to over-emphasize source-rich areas, which is a result of the summation of all sensitivities along the ray paths emitted from the common source, while obscuring areas with a low density of overlapping ray paths. On the other hand, the double-difference approach, by joining all possible pairs of measurements, does cancel out much of the shared sensitivity and, instead, illuminates wavespeed structure locally among station clusters. With a more realistic geometry on the sphere, we would expect to recover structure in the polar regions and at longitude $180^\circ$, as well as gain slight improvements for the interior structure.

Realistic phase speed input models, regional seismic arrays
-----------------------------------------------------------

Perhaps the greatest contribution of the double-difference tomographic technique as we developed it here will lie in its application to regional seismic arrays. The pairwise combination of seismic measurements should enable imaging local structures with high spatial resolution.

In Experiment V we continue the test of the synthetic model used in the global-scale Experiment IV, but this time with station located at the sites of the Transportable Array (TA) system in the United States. A view of the target model with the source-station geometry marked is shown in Fig. \[global\_test2\]a. The average station spacing is about 0.7$^\circ$. The double-difference technique is expected to resolve scale lengths comparable to the station spacing, which is much smaller than the scale of structural variations in this model. We randomly selected 300 stations from the USArray station pool for tomographic testing using membrane surface waves as in all prior experiments.

The initial model is homogeneous. The inversion result in Fig. \[global\_test2\]b was obtained from the adjoint tomography with conventional cross-correlation traveltime measurements recorded at all stations. The model resolves little more than the average wavespeed structure along the common ray path between the earthquake and the array. In contrast, the double-difference method reveals substantial details of the local structure beneath the array, and the perturbations along the direct path are also reasonable. The intrinsic sensitivity of differential measurements will make the double-difference technique ideal for the high-resolution investigation of well-instrumented areas with limited natural seismic activity, or where destructive sources are not an option.

C O N C L U S I O N S
=====================

In conventional traveltime tomography, source uncertainties and systematic errors introduce artifacts or bias into the estimated wavespeed structure. In the double-difference approach, measurements made on pairs of stations that share a similar source and similar systematic effects are differenced to partially cancel out such consequences. Both theoretically and in numerical experiments, we have shown that the ‘double-difference’ approach produces reliable misfit kernels for adjoint tomography, despite potentially incorrect source wavelets, scale factors, or timings used in synthetic modeling.

Furthermore, double-difference adjoint tomography is an efficient technique to resolve local structure with high resolution. By pairing seismograms from distinct stations with a common source, the technique yields sharp images in the area local to the stations. Those areas are often plagued by smearing artifacts in conventional tomography, especially when limited numbers of earthquakes are available nearby. This opens a wide range of promising applications in seismic array analysis, and especially in areas where the use of active sources is restricted. Fine-scale structure within the array area can be resolved even from a single earthquake.

The double-difference algorithm can be simply implemented into a conventional adjoint tomography scheme, by pairing up all regular measurements. The adjoint sources necessary for the numerical gradient calculations are obtained by summing up contributions from all pairs, which are simultaneously back-propagated to calculate the differential misfit kernels using just one adjoint simulation per iteration.

A c k n o w l e d g e m e n t s
===============================

YOY expresses her thanks to Min Chen, Paula Koelemeijer, Guust Nolet, and Youyi Ruan for discussions. This research was partially supported by NSF grants EAR-1150145 to FJS and EAR-1112906 to JT, and by a Charlotte Elizabeth Procter Fellowship from Princeton University to YOY. We thank Carl Tape and Vala Hjörleifsdóttir for sharing unpublished notes related to their published work which we cite. The Associate Editor, Jean Virieux, Carl Tape, and one other, anonymous, reviewer are thanked for their thoughtful and constructive comments, which improved our manuscript. Our computer codes are available from <https://github.com/yanhuay/seisDD>.

A P P E N D I X:D O U B L E - D I F F E R E N C EM U L T I T A P E RM E A S U R E M E N T S
===========================================================================================

In the main text we focused on measuring traveltime shifts via cross-correlation, a time-domain method that does not offer frequency resolution beyond what can be achieved using explicit filtering operations prior to analysis. Such an approach is implicitly justified for the inversion of body waves, but for strongly dispersive waves such as surface waves, we recommend making frequency-dependent measurements of phase and amplitude using the multitaper method [@Park+87a; @Laske+96; @Zhou+2004; @Tape+2009]. In this section we discuss its application to the double-difference framework. Our derivations borrow heavily from [@Hjorleifsdottir2007]. They are purposely ‘conceptual’ for transparency, and to remain sufficiently flexible for a variety of practical implementation strategies.

Multitaper measurements of differential phase and amplitude
-----------------------------------------------------------

For notational convenience we distinguish time-domain functions such as $s(t)$ or $d(t)$ from their frequency-domain counterparts $s{(\omega)}$ and $d{(\omega)}$ only through their arguments $t$ for time, and $\omega$ for angular frequency. For two stations $i$ and $j$ (always subscripted and thus never to be confused with the imaginary number $i=\sqrt{-1}$), the Fourier-transformed waveforms are assumed to be of the form [e.g., @Dziewonski+72; @Kennett2002b]

\[bli\]$$\begin{aligned}
    \label{eq:FT_si}    s_i{(\omega)}&= A_i{^{\mathrm{syn}}}{(\omega)}\,e^{-i\omega \,t_i{^{\mathrm{syn}}}{(\omega)}} 
    {\qquad\mbox{and}\qquad}s_j{(\omega)}= A_j{^{\mathrm{syn}}}{(\omega)}\,e^{-i\omega\, t_j{^{\mathrm{syn}}}{(\omega)}}
    ,\\
    \label{eq:FT_di}    d_i{(\omega)}&= A_i{^{\mathrm{obs}}}{(\omega)}\,e^{-i\omega \,t_i{^{\mathrm{obs}}}{(\omega)}}
    {\qquad\mbox{and}\qquad}d_j{(\omega)}= A_j{^{\mathrm{obs}}}{(\omega)}\,e^{-i\omega\, t_j{^{\mathrm{obs}}}{(\omega)}}
    .
  \end{aligned}$$

We define the differential frequency-dependent ‘phase’ (compare with eqs \[eq:cc\]) and ‘logarithmic amplitude’ anomalies as

\[pa\] $$\begin{aligned}
    \label{eq:def_tau_syn_ij}    {\Delta t{_{ij}}{^{\mathrm{syn}}}}{(\omega)}&= t_i{^{\mathrm{syn}}}{(\omega)}-t_j{^{\mathrm{syn}}}{(\omega)}{\qquad\mbox{and}\qquad}{\Delta\!\ln\!{A}}{_{ij}}{^{\mathrm{syn}}}{(\omega)}= \ln\left[\frac{A_i{^{\mathrm{syn}}}{(\omega)}}{A_j{^{\mathrm{syn}}}{(\omega)}}\right]
    ,
    \\ 
    \label{eq:def_tau_obs_ij}    {\Delta t{_{ij}}{^{\mathrm{obs}}}}{(\omega)}&= t_i{^{\mathrm{obs}}}{(\omega)}-t_j{^{\mathrm{obs}}}{(\omega)}{\qquad\mbox{and}\qquad}{\Delta\!\ln\!{A}}{_{ij}}{^{\mathrm{obs}}}{(\omega)}= \ln\left[\frac{A_i{^{\mathrm{obs}}}{(\omega)}}{A_j{^{\mathrm{obs}}}{(\omega)}}\right]
    .
  \end{aligned}$$

We adopt the idealized viewpoint that a trace recorded at a station $i$ may be mapped onto one made at a station $j$ [perhaps after specific phase-windowing, see @Maggi+2009 and/or after time-shifting and scaling by a prior all-frequency cross-correlation measurement] via a ‘frequency-response’ or ‘transfer’ function, ${\mathcal{T}}{_{ij}}{(\omega)}$,

\[blu\] $$\begin{aligned}
\label{eq:tf_syn} 
    s_i{(\omega)}&={{\mathcal{T}}{_{ij}}{^{\mathrm{syn}}}}{(\omega)}\,s_j{(\omega)},\\
    \label{eq:tf_obs} 
    d_i{(\omega)}&={\mathcal{T}}{_{ij}}{^{\mathrm{obs}}}{(\omega)}\,d_j{(\omega)}.
  \end{aligned}$$

Combining eqs (\[bli\])–(\[blu\]), we can rewrite the theoretical ‘transfer functions’ of our model as

\[junk\] $$\begin{aligned}
    \label{eq:tf_syn_ij}
    {{\mathcal{T}}{_{ij}}{^{\mathrm{syn}}}}{(\omega)}&=
    e^{{\Delta\!\ln\!{A}}{_{ij}}{^{\mathrm{syn}}}{(\omega)}-i\omega {\Delta t{_{ij}}{^{\mathrm{syn}}}}{(\omega)}}
    ,\\
    \label{eq:tf_obs_ij} 
    {\mathcal{T}}{_{ij}}{^{\mathrm{obs}}}{(\omega)}&=
    e^{{\Delta\!\ln\!{A}}{_{ij}}{^{\mathrm{obs}}}{(\omega)}-i\omega {\Delta t{_{ij}}{^{\mathrm{obs}}}}{(\omega)}}
    .
  \end{aligned}$$

We have hereby defined ${\Delta t{_{ij}}{^{\mathrm{syn}}}}{(\omega)}$ and ${\Delta t{_{ij}}{^{\mathrm{obs}}}}{(\omega)}$ as the frequency-dependent phase differences between $s_i(t)$ and $s_j(t)$, and $d_i(t)$ and $d_j(t)$ respectively, and ${\Delta\!\ln\!{A}}{_{ij}}{^{\mathrm{syn}}}{(\omega)}$ and ${\Delta\!\ln\!{A}}{_{ij}}{^{\mathrm{obs}}}{(\omega)}$ with their frequency-dependent amplitude differences.

The linear time-invariant model (\[blu\]) is far from perfect [for a time-and-frequency dependent viewpoint, see, e.g., @Kulesh+2005; @Holschneider+2005]: the transfer functions must be ‘estimated’. Following [@Laske+96], we begin by multiplying each of the individual time-domain records $s_i(t)$ and $d_i(t)$ by one of a small set of $K$ orthogonal data windows denoted $h^{T\Omega}_k(t)$ [here, the prolate spheroidal functions of @Slepian78], designed for a particular record length $T$ and for a desired resolution angular-frequency half-bandwidth $\Omega$ [@Mullis+91; @Simons+2015], and Fourier-transforming the result to yield the sets $$s_i^k{(\omega)}, s_j^k{(\omega)}{\qquad\mbox{and}\qquad}d_i^k{(\omega)}, d_j^k{(\omega)}{\quad\mbox{where}\quad}k=1,\ldots,\left\lfloor \frac{T\Omega}{\pi}
  \right\rfloor-1
  .$$ At frequencies separated by about $2\Omega$, any differently ‘tapered’ pairs $s^k,s^{l}$ or $d^k,d^{l}$ will be almost perfectly uncorrelated statistically [@Percival+93]. Using all of them, the transfer function in eq. (\[eq:tf\_syn\_ij\]) can be estimated in the least-squares sense [@Bendat+2000 Ch. 9] as minimizing the linear-mapping ‘noise’ power associated with eq. (\[eq:tf\_syn\]), $$\label{test}
  {T{_{ij}}{^{\mathrm{syn}}}}{(\omega)}=
  \frac{{\sum_k^K}s_i^k{(\omega)}\big[s^k_j{(\omega)}\big]^*}
       {{\sum_k^K}s^k_j{(\omega)}\big[s^k_j{(\omega)}\big]^*}=
       \operatorname*{\arg\!\min}_{{{\mathcal{T}}{_{ij}}{^{\mathrm{syn}}}}}{\sum_k^K}\left\{s_i^k{(\omega)}-{{\mathcal{T}}{_{ij}}{^{\mathrm{syn}}}}{(\omega)}s^k_j{(\omega)}\right\}^2 
       .$$ The superscripted asterisk ($^*$) denotes complex conjugation. We note here for completeness that the associated ‘coherence’ estimate is $$\label{coh}
  \gamma^2{_{ij}}{(\omega)}=
  \frac{\left|{\sum_k^K}s_i^k{(\omega)}\big[s^k_j{(\omega)}\big]^*\right|^2}
       {{\sum_k^K}s^k_i{(\omega)}\big[s^k_i{(\omega)}\big]^*
         {\sum_l^K}s^l_j{(\omega)}\big[s^l_j{(\omega)}\big]^*}
       .$$ Switching to polar coordinates we write the estimate of the transfer function in terms of its ‘gain factor’ and ‘phase angle’ contributions, $$\label{polar}
  {T{_{ij}}{^{\mathrm{syn}}}}{(\omega)}=\big|{T{_{ij}}{^{\mathrm{syn}}}}{(\omega)}\big|e^{-i\phi{_{ij}}{^{\mathrm{syn}}}(\omega)}
  ,$$ where $$\label{jobo}
  \big|{T{_{ij}}{^{\mathrm{syn}}}}{(\omega)}\big|=
  \left(\big\{{\operatorname{Im}}\big[{T{_{ij}}{^{\mathrm{syn}}}}{(\omega)}\big]\big\}^2+
  \big\{{\operatorname{Re}}\big[{T{_{ij}}{^{\mathrm{syn}}}}{(\omega)}\big]\big\}^2\right)^{1/2}
  {\qquad\mbox{and}\qquad}\phi{_{ij}}{^{\mathrm{syn}}}(\omega)=-{\mbox{atan}}\left\{
  \frac
      {{\operatorname{Im}}\big[{T{_{ij}}{^{\mathrm{syn}}}}{(\omega)}\big]}
      {{\operatorname{Re}}\big[{T{_{ij}}{^{\mathrm{syn}}}}{(\omega)}\big]}
      \right\}
      .$$ Uncertainties on the phase angle, gain factor, and coherence estimates are approximately [@Bendat+2000; @Simons+2003a] $$\label{thunc}  {\mbox{var}}\{\phi{_{ij}}{^{\mathrm{syn}}}(\omega)\}
  \approx\frac{1-\gamma^2{_{ij}}{(\omega)}}{2K{\hspace{0.1em}}\gamma^2{_{ij}}{(\omega)}}
  $$$$,\qquad
  {\mbox{var}}\big\{|{T{_{ij}}{^{\mathrm{syn}}}}|\big\}
  \approx\frac{1-\gamma^2{_{ij}}{(\omega)}}{2K{\hspace{0.1em}}\gamma^2{_{ij}}{(\omega)}}\big|{T{_{ij}}{^{\mathrm{syn}}}}\big|^2
  {\qquad\mbox{and}\qquad}{\mbox{var}}\big\{\gamma^2{_{ij}}{(\omega)}\}
  \approx
  \frac{\big[1-\gamma^2{_{ij}}{(\omega)}\big]^2}{K/2}\gamma^2{_{ij}}{(\omega)},$$ where we have silently subscribed to the established practice of confusing the estimated (on the left hand side of the expressions) from the ‘true’ quantities (on the right). Alternative estimates of the quantities in eqs (\[test\]) and (\[coh\]) could be obtained through ‘jackknifing’ [@Chave+87; @Thomson+91]. We follow [@Laske+96] in using eq. (\[test\]) to estimate the transfer function, but we do estimate its variance through the jackknife, which will be conservative [@Efron+81]. More detailed distributional considerations are provided by [@Carter87]. The utility of eq. (\[thunc\]) lies in the interpretation that highly ‘coherent’ (in the sense of eq. \[coh\]) waveforms estimated from a relatively large number of multitapered Fourier-domain ‘projections’ will be, quite generally, very precise.

Through eqs. (\[eq:tf\_syn\_ij\]) and (\[polar\]), finally, we arrive at our estimates for the frequency-dependent differential phase and amplitude measurements on the synthetic time-series $s_i(t)$ and $s_j(t)$, namely

\[stuff\] $$\begin{aligned}
    \label{eq:tau_syn_ij}    \Delta t{^{\mathrm{syn}}}{_{ij}}{(\omega)}&=
    \phi{_{ij}}{^{\mathrm{syn}}}(\omega)/\omega
    {\qquad\mbox{and}\qquad}{\Delta\!\ln\!{A}}{_{ij}}{^{\mathrm{syn}}}{(\omega)}=
    \ln\big|{{T{_{ij}}{^{\mathrm{syn}}}}{(\omega)}}\big|
    . 
  \end{aligned}$$ Exchanging the synthetic time-series for the observed ones in the above expressions, $s\leftrightarrow d$, accomplishes the desired switch in the superscripted $\mathrm{syn}\leftrightarrow\mathrm{obs}$, whence $$\begin{aligned}
    \label{eq:tau_obs_ij}    \Delta t{^{\mathrm{obs}}}{_{ij}}{(\omega)}&=
    \phi{_{ij}}{^{\mathrm{obs}}}(\omega)/\omega
    {\qquad\mbox{and}\qquad}{\Delta\!\ln\!{A}}{_{ij}}{^{\mathrm{obs}}}{(\omega)}=
    \ln\big|{T{_{ij}}{^{\mathrm{obs}}}{(\omega)}}\big|
    . 
  \end{aligned}$$

Error propagation will turn eqs (\[thunc\]) into the relevant uncertainties on the differential phase and amplitude measurements in eqs (\[stuff\]).

The logarithmic amplitude measurements in eqs (\[stuff\]) are not skew-symmetric under a change of the station indices, since eq. (\[test\]), the estimated transfer function from station $i$ to $j$, is not simply the inverse of the transfer function from station $j$ to $i$, as indeed it need not be: the modeling imperfection in mapping the seismograms is not simply invertible. Instead we find from eqs (\[test\])–(\[coh\]) that $T{_{ij}}T{_{ji}}=\gamma^2{_{ij}}$, suggesting that the required symmetry prevails only at frequencies at which the waveforms are highly coherent. A transfer-function estimate alternative to eq. (\[test\]) might thus be $T{_{ij}}/|\gamma{_{ij}}|$, as under this scenario a switch in the indices would indeed lead to ${\Delta\!\ln\!{A}}{_{ji}}$ being equal to $-{\Delta\!\ln\!{A}}{_{ij}}$ under all circumstances of mutual coherency. In the least-squares interpretation we would then be minimizing the ‘noise’ power after mapping one seismogram onto another subject to weighting by the absolute coherence to avoid mapping frequency components at which the signals are incoherent. We do not make this choice [or indeed, any other of many possible options listed by, e.g., @Knapp+76] here, but we will revisit the issue of symmetry later in this Appendix.

If the seismogram pairs had been pre-aligned and pre-scaled using conventional cross-correlation, perhaps to improve the stability of the frequency-dependent phase and amplitude expressions (\[jobo\]), these time-domain, frequency-independent, measurements would need to be added to the traveltime and amplitude expressions (\[stuff\]). If specific time windows of interest were being involved, all of the subsequent formalism would continue to hold but, naturally, the effective number of measurement pairs would increase proportionally. And of course, the double-difference multitaper method could follow any number of prior inversions designed to bring the starting model closer to the target; in this context we mention our own work on wavelet multiscale waveform-difference, envelope-difference and envelope-cross-correlation adjoint techniques [@Yuan+2014a; @Yuan+2015].

Double-difference multitaper measurements, misfit functions, and their gradients
--------------------------------------------------------------------------------

Making differential measurements between stations resolves the intrinsic ambiguity in the attribution of the phase and amplitude behavior of individual seismic records to source-related, structure-induced or instrumental factors [@Kennett2002b], and for this reason has been the classical approach for many decades [@Dziewonski+69; @Dziewonski+72]. But rather than using the differential measurements for the interrogation of subsurface structure on the interstation portion of the great circle connecting one source to one pair of receivers, we now proceed to making, and using, differential measurements between all pairs of stations into a tool for adjoint tomography. We take our cue from eq. (\[eq:cc\_DD\]) and use eqs (\[stuff\]) to define the ‘double-difference’ frequency-dependent traveltime and amplitude measurements made between stations $i$ and $j$, namely

$$\begin{aligned}
    {\Delta\Delta}t{_{ij}}{(\omega)}&={\Delta t{_{ij}}{^{\mathrm{syn}}}}{(\omega)}-{\Delta t{_{ij}}{^{\mathrm{obs}}}}{(\omega)}{\qquad\mbox{and}\qquad}{\Delta{\Delta\!\ln\!{A}}}{_{ij}}{(\omega)}={\Delta\!\ln\!{A}}{_{ij}}{^{\mathrm{syn}}}{(\omega)}-{\Delta\!\ln\!{A}}{_{ij}}{^{\mathrm{obs}}}{(\omega)}.
  \end{aligned}$$

As in eq. (\[eq:misfit\_CC\_DD\]), we incorporate those into a phase and an amplitude misfit function defined over the entire available bandwidth and for all station pairs, i.e.

\[ddmisfit\] $$\begin{aligned}
\label{eq:misfit_MTP_DD}
    \chi{_\Phi}{^{\tiny \mathrm{dd}}}&= {\frac{1}{2}}{\sum_{i}\sum_{j>i}}{\sum_\omega}W{_\Phi}{(\omega)}\big[{\Delta\Delta}t{_{ij}}{(\omega)}\big]^2,\\
    \label{eq:misfit_MTQ_DD}
    \chi{_A}{^{\tiny \mathrm{dd}}}&={\frac{1}{2}}{\sum_{i}\sum_{j>i}}{\sum_\omega}W{_A}{(\omega)}\big[{\Delta{\Delta\!\ln\!{A}}}{_{ij}}{(\omega)}\big]^2
    ,
  \end{aligned}$$

with $W{_\Phi}{(\omega)}$ and $W{_A}{(\omega)}$ appropriate weight functions. The weighted sum over all frequencies is meant to be respectful of the effective bandwidth of the estimates derived via eq. (\[test\]) from seismic records $s_i(t)$, $s_j(t)$ and $d_i(t)$, $d_j(t)$ multitapered by the Slepian window set $h^{T\Omega}_k(t)$, $k=1,\ldots,\left\lfloor T\Omega/\pi \right\rfloor-1$. The multitaper procedure implies ‘tiling’ the frequency axis over the entire bandwidth in pieces of phase and amplitude ‘information’ extracted over time portions of length $T$, centered on individual frequencies $\omega$, that are approximately uncorrelated when spaced $2\Omega$ apart [@Walden90a; @Walden+95]. The simplest form of the weight function is thus $W{_\Phi}{(\omega)}=W{_A}{(\omega)}=2\Omega$, for a discrete set of angular frequencies $\omega=\Omega,3\Omega,\ldots$ until the complete signal bandwidth is exhausted, tapering slightly at the edges of the frequency axis to smoothly contain the full bandwidth of the seismograms. Alternatively, to serve the purpose of filtering out frequencies containing incoherent energy [@Bendat+2000 Ch. 6], $W{_\Phi}{(\omega)}$ and $W{_A}{(\omega)}$ can be the inverse of the phase and amplitude variances (\[thunc\]), and, if subsampling is undesirable, they should reflect the inverse full covariances of the estimates between frequencies [@Percival+93].

As in eq. (\[eq:derivative\_misfit\_CC\_DD\]), the derivatives of the differential objective functions in eqs (\[ddmisfit\]) are

\[bliblo\] $$\begin{aligned}
    \delta \chi{_\Phi}{^{\tiny \mathrm{dd}}}&=\label{eq:derivative_MTP_DD}
    {\sum_{i}\sum_{j>i}}{\sum_\omega}W{_\Phi}{(\omega)}[{\Delta\Delta}t{_{ij}}{(\omega)}]
          \,\delta\Delta t{^{\mathrm{syn}}}{_{ij}}{(\omega)},\\
          \delta \chi{_A}{^{\tiny \mathrm{dd}}}&=\label{eq:derivative_MTQ_DD}
          {\sum_{i}\sum_{j>i}}{\sum_\omega}W{_A}{(\omega)}[{\Delta{\Delta\!\ln\!{A}}}{_{ij}}{(\omega)}]
                \,{\delta{\Delta\!\ln\!{A}}}{^{\mathrm{syn}}}{_{ij}}{(\omega)},
  \end{aligned}$$

where $\delta \Delta t{^{\mathrm{syn}}}{_{ij}}{(\omega)}$ is the perturbation of the frequency-dependent differential traveltime ${\Delta t{_{ij}}{^{\mathrm{syn}}}}{(\omega)}$ due to a model perturbation, and $\delta {\Delta\!\ln\!{A}}{^{\mathrm{syn}}}{_{ij}}{(\omega)}$ the perturbation of the frequency-dependent differential amplitude ${\Delta\!\ln\!{A}}{_{ij}}{^{\mathrm{syn}}}{(\omega)}$. It is for the terms $\delta \Delta
t{^{\mathrm{syn}}}{_{ij}}{(\omega)}$ and $\delta {\Delta\!\ln\!{A}}{^{\mathrm{syn}}}{_{ij}}{(\omega)}$ that we will construct expressions leading to a new set of adjoint sources.

Double-difference multitaper adjoint sources
--------------------------------------------

By taking perturbations to the formula in (\[eq:tf\_syn\_ij\]) in its estimated (italicized) form, we have $$\begin{aligned}
\label{eq:tf_syn_perturb}
  \delta {T{_{ij}}{^{\mathrm{syn}}}}{(\omega)}&=   {T{_{ij}}{^{\mathrm{syn}}}}{(\omega)}\big[ \delta {\Delta\!\ln\!{A}}{^{\mathrm{syn}}}{_{ij}}{(\omega)}- i\omega \delta \Delta t{^{\mathrm{syn}}}{_{ij}}{(\omega)}\big].\end{aligned}$$ Dropping the argument that captures the dependence on frequency for notational convenience, we use the product and chain rules for differentiation to write the total derivative of eq. (\[test\]) as $$\begin{aligned}
\label{eq:tf_syn_mt_perturb}
  \delta {T{_{ij}}{^{\mathrm{syn}}}}{}&=\frac{1}{S{_{jj}}}{\sum_k^K}\, [s_j^k{}]^* \,\delta s_i^k{}-{T{_{ij}}{^{\mathrm{syn}}}}{}[s_j^k{}]^* \,\delta s_j^k{}+\left( s_i^k{}-{T{_{ij}}{^{\mathrm{syn}}}}{}s_j^k{}\right)
  \,\delta [s_j^k{}]^*
  .\end{aligned}$$ We introduced a symbol for the multitapered (cross-) power-spectral density estimates $$S{_{ij}}={\sum_k^K}s^k_i{}[s^k_j{}]^*
  {\quad\mbox{with}\quad}{T{_{ij}}{^{\mathrm{syn}}}}={\frac{\displaystyle{S{_{ij}}}}{\displaystyle{S{_{jj}}}}}
  {\qquad\mbox{and}\qquad}S{_{ij}}^*=S{_{ji}},$$ which of course upholds the relation eq. (\[test\]). Combining eqs. (\[eq:tf\_syn\_perturb\])–(\[eq:tf\_syn\_mt\_perturb\]) and using eq (\[test\]) again, we then obtain the perturbations $$\begin{aligned}
\label{eq:dt_perturb}
  \delta \Delta t{^{\mathrm{syn}}}{_{ij}}{(\omega)}&=\frac{i}{2\omega}
  {\sum_k^K}\frac{[s_j^k{}]^*}{S{_{ij}}} \,\delta s_i^k{}-\frac{s_j^k{}}{S{_{ij}}^*} \,\delta [s_i^k{}]^* 
  -\frac{[s_i^k{}]^*}{S{_{ij}}^*} \,\delta s_j^k{}+\frac{s_i^k{}}{S{_{ij}}}\,\delta [s_j^k{}]^*\end{aligned}$$ through straightforward manipulation, and, likewise, $$\begin{aligned}
\label{eq:dA_perturb}
  {\delta{\Delta\!\ln\!{A}}}{^{\mathrm{syn}}}{_{ij}}{(\omega)}&={\frac{1}{2}}{\sum_k^K}\frac{[s_j^k{}]^*}{S{_{ij}}} \,\delta s_i^k{}+\frac{s_j^k{}}{S{_{ij}}^*} \,\delta [s_i^k{}]^*
  +\left(\frac{[s_i^k{}]^*}{S{_{ij}}^*}-2\frac{[s_j^k{}]^*}{S{_{jj}}}\right) \delta s_j^k{}+\left(\frac{s_i^k{}}{S{_{ij}}}-2\frac{s_j^k{}}{S{_{jj}}}\right) \delta [s_j^k{}]^*
  .\end{aligned}$$ We shall associate the frequency-dependent partials that appear in eqs (\[eq:dt\_perturb\])–(\[eq:dA\_perturb\]) with the symbols

$$\begin{aligned}
\label{ps}
    p_{i|j}^k{(\omega)}&= \frac{i}{2\omega}\frac{[s_j^k{}]^*}{S{_{ij}}}
    {\qquad\mbox{and}\qquad}p_{j|i}^k{(\omega)}= \frac{-i}{2\omega}\frac{[s_i^k{}]^*}{S{_{ij}}^*},
    \\\label{as}
    a_{i|j}^k{(\omega)}&={\frac{1}{2}}\frac{[s_j^k{}]^*}{S{_{ij}}}
    {\qquad\mbox{and}\qquad}a_{j|i}^k{(\omega)}={\frac{1}{2}}\frac{[s_i^k{}]^*}{S{_{ij}}^*}-\frac{[s_j^k{}]^*}{S{_{jj}}},
  \end{aligned}$$

and their conjugates. Again, we remark on the lack of station symmetry in the amplitude terms (\[as\]). However, it is noteworthy that, via the Fourier identity ${\partial}_t\leftrightarrow i\omega$, the terms in the expressions (\[ps\]) are identical (up to the sign) to the (first) terms in eq. (\[as\]) if in the latter, the time-domain seismograms are substituted for their time-derivatives. Consequently, we may recognize eqs (\[eq:dt\_perturb\])–(\[eq:dA\_perturb\]) as the (differential) frequency-dependent (narrow-band) multitaper generalization of the traveltime expressions (\[dti\]) and (\[dtss\]), familiar from [@Dahlen+2000], and of their amplitude-anomaly counterparts, see [@Dahlen+2002]. With these conventions in place, the derivatives of the multitaper objective functions in eq. (\[bliblo\]) are now given by

\[blonki\] $$\begin{aligned}
\label{eq:blonki}
    \delta\chi{^{\tiny \mathrm{dd}}}{_\Phi}=&
    {\sum_{i}\sum_{j>i}}{\sum_k^K}{\sum_\omega}W{_\Phi}{(\omega)}[{\Delta\Delta}t{_{ij}}{(\omega)}]
    \left\{
    p_{i|j}^k{(\omega)}\,\delta s_i^k{(\omega)}+[p_{i|j}^k{(\omega)}]^*  \,\delta [s_i^k{(\omega)}]^*
    +p_{j|i}^k{(\omega)}\,\delta s_j^k{(\omega)}+[p_{j|i}^k{(\omega)}]^* \,\delta [s_j^k{(\omega)}]^*
    \right\}
    ,\\
    \delta\chi{^{\tiny \mathrm{dd}}}{_A}=&
    {\sum_{i}\sum_{j>i}}{\sum_k^K}{\sum_\omega}W{_A}{(\omega)}[{\Delta{\Delta\!\ln\!{A}}}{_{ij}}{(\omega)}] 
    \left\{
    a_{i|j}^k{(\omega)}\,\delta s_i^k{(\omega)}+[a_{i|
        j}^k{(\omega)}]^* \,\delta [s_i^k{(\omega)}]^*
    +a_{j|i}^k{(\omega)}\,\delta s_j^k{(\omega)}+[a_{j|i}^k{(\omega)}]^* \,\delta [s_j^k{(\omega)}]^*
    \right\}
    .
  \end{aligned}$$

We introduce the time-domain functions, inverse-Fourier transformed from their frequency-domain product counterparts,

$$\begin{aligned}
    \mathscr{S}_{i|j}{(t)}&=
    {\sum_k^K}h^{T\Omega}_k{(t)}\,\mathcal{F}^{-1}\!\left\{W{_\Phi}[{\Delta\Delta}t{_{ij}}]
    \,p_{i|j}^k\right\}\!{(t)},
    {\qquad\mbox{and}\qquad}\mathscr{S}_{j|i}{(t)}=
    {\sum_k^K}h^{T\Omega}_k{(t)}\,\mathcal{F}^{-1}\!\left\{W{_\Phi}[{\Delta\Delta}t{_{ij}}]
    \,p_{j|i}^k\right\}\!{(t)},
    \\
    \mathscr{A}_{i|j}{(t)}&=
    {\sum_k^K}h^{T\Omega}_k{(t)}\,\mathcal{F}^{-1}\!\left\{W{_A}[{\Delta{\Delta\!\ln\!{A}}}{_{ij}}]
    \,a_{i|j}^k\right\}\!{(t)}{\qquad\mbox{and}\qquad}\mathscr{A}_{j|i}{(t)}=
    {\sum_k^K}h^{T\Omega}_k{(t)}\,\mathcal{F}^{-1}\!\left\{W{_A}[{\Delta{\Delta\!\ln\!{A}}}{_{ij}}]
    \,a_{j|i}^k\right\}\!{(t)},
  \end{aligned}$$

with which we are able to rewrite eqs (\[blonki\]) in the time domain, via Plancherel’s relation as applicable over the record section of length $T$. We have simultaneously extricated the tapers from their seismograms, $\delta s_i^k{(t)}=\delta
[s_i^k{(t)}]^*=h^{T\Omega}_k{(t)}\,\delta s_i{(t)}$, and used the real-valuedness of the time-series to arrive at the recognizable (see eq. \[dchi\]) form

$$\begin{aligned}
\label{eq:blonkt}
    \delta\chi{^{\tiny \mathrm{dd}}}{_\Phi}=&
    {\int_0^T}\left\{
    \sum_{i} \left[4\pi
      \sum_{j>i} \mathscr{S}_{i|j}{(t)}\right]\delta s_i{(t)}+\sum_{j} \left[4\pi
      \sum_{i<j} \mathscr{S}_{j|i}{(t)}\right]\delta s_j{(t)}\right\}{\hspace{0.1em}\mathrm{d}t},
    \\
    \delta\chi{^{\tiny \mathrm{dd}}}{_A}=&
    {\int_0^T}\left\{
    \sum_{i} \left[4\pi
      \sum_{j>i} \mathscr{A}_{i|j}{(t)}\right]\delta a_i{(t)}+\sum_{j} \left[4\pi
      \sum_{i<j} \mathscr{A}_{j|i}{(t)}\right]\delta a_j{(t)}\right\}{\hspace{0.1em}\mathrm{d}t}.
  \end{aligned}$$

The multitaper phase and amplitude adjoint sources are thus $$\begin{aligned}
\label{yolo}
  f_{\Phi,i}^{\dagger}(t)&=
  4\pi \sum_{j>i} \mathscr{S}_{i|j}{(T-t)}\,\delta ({\mathbf{x}}-{\mathbf{x}}_i)
  {\qquad\mbox{and}\qquad}f_{\Phi,j}^{\dagger}(t)=
  4\pi \sum_{i<j} \mathscr{S}_{j|i}{(T-t)}\,\delta ({\mathbf{x}}-{\mathbf{x}}_j)
  ,
  \\\label{polo}
  f_{A,i}^{\dagger}(t)&=
  4\pi \sum_{j>i} \mathscr{A}_{i|j}{(T-t)}\,\delta ({\mathbf{x}}-{\mathbf{x}}_i)
  {\qquad\mbox{and}\qquad}f_{A,j}^{\dagger}(t)=
  4\pi \sum_{i<j} \mathscr{A}_{j|i}{(T-t)}\,\delta ({\mathbf{x}}-{\mathbf{x}}_j)
  .\end{aligned}$$ It is gratifying to discover that eq. (\[yolo\]) is indeed the multitaper generalization of eqs (\[eq:adj\_CC\_DD\]), which themselves generalized eq. (\[eq:adj\_CC\]). For its part, eq. (\[polo\]), ultimately, is the double-difference multitaper generalization of the amplitude-tomography adjoint source derived by, e.g., [@Tromp+2005], but which we have not given any further consideration in the main text.

Finally, we return to the issue of symmetry in the differential amplitude estimation, which can be restored by replacing all separate summations of the type $\sum_i\sum_{j>i}$ and $\sum_j\sum_{i<j}$ that appear in the expressions above by a common $\sum_{ij}$.

With these expressions we are finally able to define the multitaper misfit sensitivity kernels along the same lines as eqs (\[kp\]) and (\[kp2\]), one for the phase and one for the amplitude measurements, both in terms of the velocity perturbation — we are (yet) not concerned with variations in intrinsic attenuation [e.g., @Tian+2009; @Tian+2011; @Zhu+2013].

Numerical experiment
--------------------

Fig. \[fig:Sin\_model\]a shows a simple two-dimensional vertical *S*-wave speed model with an anomalous layer that divides the model into three segments: a top portion with an *S*-wave speed of 800 m/s, a middle with a speed of 1000 m/s, and a bottom with a speed of 800 m/s. The initial model, shown in Fig. \[fig:Sin\_model\]b, possesses the same tripartite geometry but the wavespeeds have different values: at the top 900 m/s, in the middle 800 m/s, and at the bottom 1000 m/s.

Synthetics were computed used one point-source Ricker wavelet with 40 Hz dominant frequency located at 0.5 m depth, at a horizontal distance of 100 m from the left edge of the model. Two receivers were placed at the surface in the positions shown in Fig. \[fig:Sin\_model\].

Fig. \[fig:Sin\_data\] shows the observed and synthetic records (a–b), the adjoint sources (c–d), and the double-difference frequency-dependent multitaper phase measurement (e). Fig. \[fig:Sin\_kernel\] compares the misfit sensitivity kernel of the double-difference time-domain cross-correlation measurement to the misfit sensitivity kernel of the double-difference multitaper phase measurement, which shows some additional structure.

[![image](Fig//Sinmodel_test1_adj){width="48.00000%"}]{}
