Validating Forward Modeling and Inversions of Helioseismic Holography Measurements

Helioseismology has been a useful tool for studying the subsurface properties of the Sun. It is generally divided into global (e.g., Christensen-Dalsgaard 2003; Howe 2009) and local (e.g., Gizon et al. 2010; Braun 2015) applications. The former include inferences of the radially symmetric structure of the Sun and the latitudinal and depth dependence of its internal rotation, while the latter include studies of the relatively small-scale structure and flows below sunspots, active regions, and supergranulation, as well as larger-scale convection and meridional circulation. Giving us the means to indirectly image the Sun's interior, helioseismology is of great importance to the study of solar structure and subsurface dynamics.

The forward problem in helioseismology is to determine the relationship between some quantity, measurable at the solar surface, and an unobserved subsurface feature that we would like to study. For local applications like those considered in this work, the former is helioseismic wave travel-time measurements, and the latter is vector plasma flows. The two are related through linear integral equations of the form

$$\begin{eqnarray}&&\delta {\tau }_{a}({\boldsymbol{r}})={\int }_{\odot }\displaystyle \sum _{\beta }{K}_{\beta }^{a}({\boldsymbol{r}}^{\prime} -{\boldsymbol{r}};z){v}_{\beta }({\boldsymbol{r}}^{\prime} ,z){d}^{2}{\boldsymbol{r}}^{\prime} \,{dz}+{n}^{a}({\boldsymbol{r}}),\end{eqnarray} \tag{ 1 }$$

where, for each travel-time measurement δτ_a, ${K}_{\beta }^{a}$ is the three β ∈ {x, y, z} components of a set of vector-valued functions called sensitivity kernels and n^a is the noise in these measurements. Here, ${\boldsymbol{r}}=(x,y)$ is the horizontal position and z is the height (noting that z = 0 at the surface and z < 0 inside the Sun). Given these quantities, the inverse problem is to solve for the subsurface flow v_β as accurately as possible through a series of matrix inversions.

In the past decade, the development and availability of realistic artificial data, obtained from numerical wave-propagation computations, have allowed the validation and testing of helioseismic procedures. Relevant numerical simulations include those computed under hydrostatic or magnetohydrostatic conditions (e.g., Hanasoge et al. 2007; Parchevsky & Kosovichev 2007; Hartlep et al. 2008; Felipe et al. 2010) as well as fully compressible hydrodynamic or magnetohydrodynamic computations (e.g., Rempel et al. 2009; Stein et al. 2009; Rempel 2015; Stein & Nordlund 2012). Simulations give us the opportunity to test kernels and inversion procedures on data whose flow structure is known a priori and whose properties resemble the real Sun as closely as possible. This kind of validation has been performed for both time–distance (e.g., Zhao et al. 2007; Švanda et al. 2011; DeGrave et al. 2014a, 2014b; Parchevsky et al. 2014) and helioseismic holography (HH) methods (e.g., Braun et al. 2007; Birch et al. 2009; Braun et al. 2012; Dombroski et al. 2013; Braun 2014). Validation tests of this kind are an important and necessary step in the helioseismic analysis of solar subsurface structure, as we are interested in recovering information from regions of the Sun that cannot be directly observed.

In this work, we test the performance of HH and kernels through a series of forward and inverse-modeling comparisons, employing a realistic numerical simulation of the quiet Sun. The forward and inverse-modeling tests we conduct employ a set of kernels computed for travel times measured using HH carried out in lateral-vantage geometry (Lindsey & Braun 2004). Until now, tests of lateral-vantage HH have only been carried out through comparisons of measured and forward-modeled travel times (Braun et al. 2007). This work presents the first comprehensive test of inverse modeling for flows using lateral-vantage HH.

The layout of the paper proceeds as follows: In Section 2, we describe the quiet-Sun simulation data used in this work. The holography travel-time measurement procedure is discussed in Section 3, and the forward modeling is detailed in Section 4. The inversion method and results are described in Sections 5 and 6, and concluding remarks are given in Section 7.

The simulation we employ in this work represents quiet-Sun convection with a small-scale dynamo and has been described in detail in Rempel (2014). A time series of artificial data from this simulation has been previously employed in validating the inversion procedures used in this work (DeGrave et al. 2014a), using time–distance measurements and sensitivity kernels. The simulation was computed with horizontal and vertical resolutions of 64 km and 32 km, respectively. Convective motions excite surface gravity and acoustic waves that propagate throughout the domain that spans 98.3 × 98.3 Mm horizontally and 18.4 Mm vertically. A cut through the (15 hr) time-averaged v_x flow component of the simulation is shown in Figure 1.

Zoom In Zoom Out Reset image size

Doppler velocity time series, assuming a vertical line-of-sight, extracted at an optical depth of 0.01 and having a cadence of 45 s, are publicly available.³ These artificial Dopplergrams have been interpolated from the original simulation onto a coarser grid with a horizontal spacing of 384 km. For this work we utilize the first 15 hr of the 30 hr simulation run.

HH is a method that computationally extrapolates the surface acoustic field from a selected pupil into the solar interior (Lindsey & Braun 1997) in order to estimate the complex amplitudes of the waves propagating into or out of a focus point at a chosen depth within the solar interior. These amplitudes are called the acoustic ingression and acoustic egression, respectively. Lateral-vantage holography (e.g., Lindsey & Braun 2004) is analogous to deep-focus methods in time–distance helioseismology and common-depth-point reflection terrestrial seismology. Figure 2 illustrates the pupil geometry used. The annulus is defined by rays propagating through the focus and inclined up to ±45° from the direction parallel to the surface. The practical aspects of the methodology have been described in detail elsewhere (e.g., Braun et al. 2007; Braun & Birch 2008; Braun 2014). Most of the prior applications of lateral-vantage HH make use of the northward-minus-southward (NS) and westward-minus-eastward (WE) travel-time differences δτ_ns and δτ_we as derived from cross-covariances between the egression and ingression as assessed using opposite quadrants (Figure 2(b)). Here, an additional pair of cross-covariances are obtained using inner and outer portions of the complete annulus (Figure 2(c)). The radius ρ_h that separates the two subannuli is defined by the ray path of a wave propagating horizontally through the focus. Cross-covariances between egressions and ingressions assessed in these subannuli are used to determine an outward-minus-inward (OI) propagation travel-time difference δτ_oi. As described elsewhere (Braun 2014), Gaussian phase-speed filters are also used, with a width δw and a peak at w_o corresponding to the phase speed of the aforementioned horizontally propagating wave. The minimum and maximum radii of the annulus, ρ_min and ρ_max, horizontal-ray radius ρ_h, and phase-speed filter parameters for the measurements used here are given in Table 1.

Zoom In Zoom Out Reset image size

To compute kernels for deep-focusing travel times measured from the synthetic observations, we first construct a model power spectrum by fitting the power spectrum of the synthetic data. We fit for a source function (this determines the mode amplitude), damping rate, and deviation of the frequency from Model S (Christensen-Dalsgaard et al. 1996) frequency as functions of horizontal wavenumber and radial order. The details are described in Appendix A. We then use the Born approximation to compute the sensitivity of the lateral-vantage travel-time differences to small-amplitude steady flows. The calculation is based on the approach of Birch & Gizon (2007), though extended to include the holography Greens functions and the pupil functions. The product of the Greens function and the pupils appears in the calculation in exactly the same way as a non-axisymmetric complex-valued data-analysis filter. In this calculation, we use the source function and damping rates obtained from the fit to the power spectrum. We, however, do not include in the calculations the effect of the changes in the mode frequencies but instead use the normal-mode frequencies and eigenfunctions of Model S. We estimate that this approximation introduces an error of about 2% in the kernels. For the applications in this work, this error is significantly smaller than the noise in the travel-time measurements (see Section 4.2). Appendix B shows the details of the calculation.

4.1. Model for the Noise Covariance

4.2. Forward Comparisons

Figure 3 shows comparisons of travel-time difference maps measured from the simulation using HH with maps predicted from the sensitivity functions convolved with the true time-averaged flows present in the simulation. Correlation statistics comparing each pair of measured and forward-modeled maps are shown in Table 2 for the WE and OI travel-time differences. The table includes the rms of the difference between maps (rms error) and the slope of the least-squares linear fit between measured and modeled values (the fit assumes there are no errors in the modeled values). Good agreement is found between measured and forward-modeled travel-time maps, and slope values are close to unity for all focus depths.

Zoom In Zoom Out Reset image size

Table 2.  Lateral-vantage Correlation Statistics Comparing Measured and Forward-modeled Travel-time Differences

Travel-time	Focus Depth	rms Error	Slope
Difference	(Mm)	(s)
δτwe	0.77	7.3	1.02
δτwe	1.53	5.6	0.99
δτwe	2.30	4.7	1.03
δτwe	2.99	4.3	1.00
δτwe	3.97	3.8	0.98
δτwe	5.01	3.8	0.98
δτwe	5.99	3.8	0.98
δτwe	6.96	3.8	0.98
δτwe	8.35	4.0	0.94
δτoi	0.77	3.7	1.01
δτoi	1.53	3.3	0.98
δτoi	2.30	3.3	1.01
δτoi	2.99	3.2	1.00
δτoi	3.97	3.2	1.03
δτoi	5.01	3.4	1.06
δτoi	5.99	3.2	1.04
δτoi	6.96	3.2	1.01
δτoi	8.35	3.4	0.99

Download table as:  ASCIITypeset image

The rms errors have values between 4 and 7 s for the WE measurements and between 3 and 4 s for the OI measurements. For context, the forward-modeled maps exhibit travel-time differences with peak values ranging from 60 s for the shallowest focus depth to about 15 s for the deepest measurements. We note that a 2% error in the kernels (as discussed above) would produce travel-time errors that are, for the most part, considerably less than 1 s. This is significantly smaller than the rms errors.

To recover flows from the simulation, we employ the subtractive optimally localized averaging (SOLA) method (Pijpers & Thompson 1992). The goal of SOLA is to find a set of two-dimensional inversion weights (see Švanda et al. 2011; Jackiewicz et al. 2012) that, when spatially convolved with the travel-time measurements, will give a smoothed estimate of flow component α = {x, y, z} at some target depth z₀ within the simulation domain:

$$\begin{eqnarray}&&{v}_{\alpha }^{\mathrm{inv}}({\boldsymbol{r}};{z}_{0})=\displaystyle \sum _{i}\displaystyle \sum _{a=1}^{M}{w}_{\alpha }^{a}({{\boldsymbol{r}}}_{i}-{\boldsymbol{r}};{z}_{0})\delta {\tau }_{a}({{\boldsymbol{r}}}_{i}),\end{eqnarray} \tag{ 2 }$$

where δτ_a is the set of M travel-time measurements and w^a_α is their respective weights. The sum over i is over all horizontal positions. When a set of weights has been computed, they are linearly combined with the sensitivity kernels to produce a so-called averaging kernel

$$\begin{eqnarray}&&{{ \mathcal K }}_{\alpha }^{\beta }({\boldsymbol{r}},z;{z}_{0})=\displaystyle \sum _{i}\displaystyle \sum _{a=1}^{M}{w}_{\alpha }^{a}({{\boldsymbol{r}}}_{i}-{\boldsymbol{r}};{z}_{0}){K}_{\beta }^{a}({\boldsymbol{r}}-{{\boldsymbol{r}}}_{i};z)\end{eqnarray} \tag{ 3 }$$

that effectively gives the spatial resolution of the inversion. Here, β = {x, y, z} is the three kernel components, and the sum over i represents a horizontal convolution. Ideally, the weights will be such that, for α = β, the resulting averaging kernel is well localized in three-dimensional space, closely matching a predefined (typically Gaussian) target function, T. For $\beta \ne \alpha$ the averaging kernel will ideally be small—these off-diagonal components are responsible for cross-talk (Jackiewicz et al. 2012).

For each inversion target depth z₀ and target flow direction α, we search for a set of weights that minimizes the cost function

$$\begin{eqnarray}&&{ \mathcal X }={\int }_{\odot }{[{{ \mathcal K }}_{\alpha }^{\alpha }({\boldsymbol{r}},z;{z}_{0})-T({\boldsymbol{r}},z;{z}_{0})]}^{2}\,{d}^{2}{\boldsymbol{r}}\,{dz}\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&+\ \nu \displaystyle \sum _{\beta \ne \alpha }{\int }_{\odot }{[{{ \mathcal K }}_{\alpha }^{\beta }({\boldsymbol{r}},z;{z}_{0})]}^{2}{d}^{2}{\boldsymbol{r}}\,{dz}\end{eqnarray} \tag{ 4b }$$

$$\begin{eqnarray}&&+\ \epsilon \displaystyle \sum _{a,i}{[{w}_{\alpha }^{a}({{\boldsymbol{r}}}_{i};{z}_{0})]}^{2}\end{eqnarray} \tag{ 4c }$$

$$\begin{eqnarray}&&+\ \mu \displaystyle \sum _{a,b,i,j}{w}_{\alpha }^{a}({{\boldsymbol{r}}}_{i};{z}_{0}){{\rm{\Lambda }}}_{{ab}}{w}_{\alpha }^{b}({{\boldsymbol{r}}}_{j};{z}_{0}),\end{eqnarray} \tag{ 4d }$$

where the various terms in Equation (4) represent the misfit between the averaging kernel and target function (4a); the extent to which the non-inverted-for flow components contribute to the recovered velocities, referred to as cross-talk (4b); an ad hoc term quantifying the localization of the inversion weights, referred to as the weight spread (4c); and the level of random noise in the solution (4d). Λ_ab is the noise covariance. These terms can be controlled to some degree by varying regularization parameters ν, , and μ. In practice, we perform inversions and compute their respective (4a–4d) quantities for many combinations of regularization values for a given target depth. An optimal combination of parameters is then selected from these. Example solution grids for a typical horizontal flow inversion are shown in Figure 4.

Zoom In Zoom Out Reset image size

When performing inversions for the horizontal flow components (v_x, v_y) in this work, priority is first placed on choosing an acceptable noise level of roughly 30 m s⁻¹ for each target depth. This noise level seems reasonable given the fact that flows in the upper 5 Mm of the simulation domain are of the order of 300 m s⁻¹. Of the various combinations of regularization parameters that yield the selected noise level, we choose the solution that provides a reasonable trade-off between the misfit and weight spread. This is done through a typical L-curve analysis. For horizontal flow inversions, we set ν = 0, and therefore do not regularize the amount of cross-talk, as its effects are small compared to the large-amplitude flows that we are inverting for.

Inversions for the vertical flow component, v_z, also place priority on first selecting an acceptable level of noise. The chosen noise level must be much lower than that of the horizontal flow inversions, as the amplitude of the vertical flows in the top layers of the simulation domain are roughly an order of magnitude weaker than the horizontal flows (≈15 m s⁻¹ on supergranule scales). We therefore choose a noise level of 5 m s⁻¹ in order to retain the possibility of recovering the weak vertical flow signal while balancing the misfit between averaging kernel and target function.

Unlike the horizontal flow inversions, we also now fix ν ≈ 100, in addition to varying the misfit and weight spread regularization parameters. Constraining cross-talk is necessary for vertical flow inversions; failing to do so can give rise to cross-talk that has the same amplitude or greater than the weak vertical flows for which we are inverting. This often leads to recovered flows that are strongly anticorrelated with the true flows (e.g., Zhao et al. 2007). Of the solutions that yield the specified noise level, we again choose the one that provides a reasonable trade-off between the misfit and weight spread.

The code used to perform inversions in this work has been validated previously with time–distance measurements applied to realistic simulations (DeGrave et al. 2014a, 2014b). To assess the performance of the inversions, we compare the recovered flows, ${v}_{\alpha }^{\mathrm{inv}}$, with the target flows, which represent the true flows ${v}_{\alpha }^{\mathrm{sim}}$ smoothed to the expected resolution of the inversion flow maps via convolution with the inversion target function:

$$\begin{eqnarray}&&{v}_{\alpha }^{\mathrm{tgt}}({\boldsymbol{r}},{z}_{0})={\int }_{\odot }T({\boldsymbol{r}}^{\prime} -{\boldsymbol{r}};z,{z}_{0}){v}_{\alpha }^{\mathrm{sim}}({\boldsymbol{r}}^{\prime} ,z)\,{d}^{2}{\boldsymbol{r}}^{\prime} \,{dz}.\end{eqnarray} \tag{ 5 }$$

The target solution represents the best we can hope to achieve in any particular inversion. For the vertical flow inversion, we also examine the cross-talk components,

$$\begin{eqnarray}&&{v}_{\alpha }^{\beta }({\boldsymbol{r}};{z}_{0})={\int }_{\odot }{{ \mathcal K }}_{\alpha }^{\beta }({\boldsymbol{r}}^{\prime} -{\boldsymbol{r}},z;{z}_{0}){v}_{\beta }^{\mathrm{sim}}({\boldsymbol{r}}^{\prime} ,z)\,{d}^{2}{\boldsymbol{r}}^{\prime} \,{dz},\end{eqnarray} \tag{ 6 }$$

in addition to ${v}_{\alpha }^{\mathrm{inv}}$. We note that ${v}_{\alpha }^{\mathrm{inv}}\equiv {v}_{\alpha }^{\alpha }$ when the level of noise is negligible.

6.1. Horizontal Flow Inversions

Inversions for the horizontal flow components (v_x, v_y) were carried out at depths of 1, 3, and 5 Mm below the surface of the simulation domain. Figure 5 shows maps of the resulting lateral-vantage inversions (top row), along with the smoothed simulation flows, ${v}_{x}^{\mathrm{tgt}}$, at each of these depths (bottom row). The two-dimensional Pearson correlation values between the inversion and target simulation flows are given in the lower left-hand corner of each panel. The noise level for each of the inversions is approximately 30 m s⁻¹ for all target depths. The horizontal resolution of each inversion (i.e., the horizontal FWHM of the target function) was 10 Mm for all depths. The vertical resolution of the inversions is 1.4 Mm for the target depth of 1 Mm, and 2 Mm for the deeper target depths. We find that the inversions are able to recover the simulation flows quite well, particularly in the upper 3 Mm of the domain, and correlation values are generally high here (>0.8). However, at the 5 Mm depth, the quality of our inversion has deteriorated significantly, and we are not able to accurately recover flows there. Flow amplitudes are well reproduced at the two shallowest depths, but are underestimated at the 5 Mm depth by a factor of roughly 2.4 in terms of (v_x, v_y) rms values.

Zoom In Zoom Out Reset image size

One-dimensional cuts through the x-component of the inversion averaging kernels and target functions are shown in Figure 6. These figures show the depths targeted in the inversions along with the depths that have actually been sampled. The averaging kernels show some undesirable characteristics, most notably the large misfit observed at z > −1 Mm for the 1 Mm depth (left panel), the strong near-surface contribution at z > −1.5 Mm for the 3 Mm depth (middle panel), and the negative near-surface lobe for the 5 Mm depth (right panel). The unwanted near-surface sensitivity, particularly at the 3 and 5 Mm depths, appears to have little effect on the recovered flows shown in Figure 5. This is likely, at least to some degree, due to the fact that the simulation flow field does not vary rapidly with depth. For example, the large-scale convective features in the simulation are very extended in depth with flows that do not show any sudden reversal in sign (e.g., outflow to inflow) or large changes in amplitude in the near-surface layers. If that were not the case, or if the true flow structure was not known a priori, these misfit issues could make it difficult to properly interpret the results, particularly for inversions in the deeper layers where flows are more difficult to retrieve.

Zoom In Zoom Out Reset image size

6.2. Vertical Flow Inversions

An inversion was also carried out to retrieve the vertical flow component, v_z, at a depth of 1 Mm. Following Švanda (2013), this inversion employs only kernels computed for the OI measurements. We neglect the WE and NS measurements as they are relatively insensitive to vertical flows and can actually adversely effect the recovered flows by contributing additional noise. The recovered flow map is shown in Figure 7. The noise level for these inversion is 5 m s⁻¹, and the spatial resolution is the same as for the 1 Mm v_x case. The two-dimensional Pearson correlation value between the inversion and target simulation flows is 0.56. Figure 8 shows cuts in depth through the averaging kernel components for the inversion. It appears that the inversion is able to minimize the cross-talk well, with the ${{ \mathcal K }}_{z}^{x}$ and ${{ \mathcal K }}_{z}^{y}$ cross-terms accounting for about 10% or less of the maximum kernel amplitude. The averaging kernel is strongly peaked near the surface of the simulation domain, but shows some extended sensitivity over a range of depths not actually targeted in the inversion. The ${{ \mathcal K }}_{z}^{z}$ component also exhibits broad negative side-lobes that we are not able to minimize effectively in the inversion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 shows the v_z inversion flow map (leftmost column) along with the contributions from the individual terms in Equation (6) (i.e., a convolution of the Figure 8 averaging kernel components with the true simulation flows). We see again that the cross-talk minimization of ${v}_{z}^{x}$ and ${v}_{z}^{y}$ is reasonably effective, but somewhat less so than Figure 8 indicates due to the strong horizontal flows with which ${{ \mathcal K }}_{z}^{x}$ and ${{ \mathcal K }}_{z}^{y}$ have been convolved. These cross-terms have rms amplitudes that are about 50% that of ${v}_{z}^{z}$. The sum of columns 2–4 is shown in column 5 and closely resembles the target simulation flows. This suggests that the cross-talk is not detrimental to the inversions in terms of being able to reproduce the correct flow structure of the simulation.

Zoom In Zoom Out Reset image size

Column 6 shows the residual of columns 1 and 5 and represents the contribution of noise to the recovered flows (Jackiewicz et al. 2012). Comparing these three columns shows that the vertical flow inversions are clearly dominated by this noise. Though we fail to adequately recover the vertical flow component here, it is important to note that these results are not dissimilar from the vertical flow inversions by DeGrave et al. (2014a) using kernels computed under the Birch & Gizon (2007) prescription.

We have introduced and successfully tested a set of sensitivity kernels for use in local helioseismology by employing them in a series of forward and inverse-modeling comparisons using HH measurements. Measured travel times computed in the lateral vantage compared favorably with forward-modeled ones predicted by the kernels, both in terms of spatial distribution and amplitude. Inversions for the horizontal flow components (v_x, v_y) employing the kernels were successful in recovering the simulation flow field from the upper 3 Mm of the domain, and flow amplitudes agreed well with those of the target flows. However, inversions carried out at a depth of 5 Mm were less successful in reproducing the flows than in the near-surface layers, and flow amplitudes were underestimated there by a factor of roughly 2.4. It is important to note, though, that the inability of the inversions to recover flows at this depth is more a consequence of noise issues rather than a problem with the kernels themselves. A near-surface inversion for the vertical flow component failed to adequately retrieve the simulation v_z flows. Though recovered flows correlated reasonably well with the target simulation flows, amplitudes were not well reproduced, and the inversion was dominated by noise.

This work represents the first comprehensive test of the ability of helioseismic holography, as employed in the lateral-vantage (deep-focus) configuration, to infer subsurface flows on spatial scales on the order of, and smaller than, supergranules. As such, it extends and confirms the general findings of a prior validation study (Braun et al. 2007) carried out for lateral-vantage HH, but using only forward-model comparisons. Specifically, Braun et al. (2007) concluded from a consideration of signal to noise (and excluding inversion-related issues) that supergranule-sized flows are undetectable below about 5 Mm using data spanning less than the lifetime of a typical supergranule. The results found here also complement the inverse-modeling validation study performed using surface-focus holography (Dombroski et al. 2013), which employed regularized-least-squares (RLS) inversions of travel times measured from simulations of an idealized supergranule-like flow.

Our results are also consistent with findings from validation studies of time–distance helioseismology (e.g., Zhao et al. 2007; Švanda et al. 2011; DeGrave et al. 2014a). It is now readily apparent that methods that explicitly include the minimization of cross-talk effects such as those presented here (and, e.g., Švanda et al. 2011; DeGrave et al. 2014a) offer distinct improvements in the determination of vertical flows over methods that do not, such as the RLS inversion of Dombroski et al. (2013), or inversions based on ray theory (Zhao et al. 2007).

The dominance of realization noise for target depths below a few Mm has lead to statistical approaches to inferring deeper flows. A notable example is the averaging of measurements made with respect to thousands of supergranules (e.g., Švanda 2012; Duvall & Hanasoge 2013; Duvall et al. 2014). Even so, inverse modeling of these and other data apparently remains challenging (e.g., Švanda 2015; Bhattacharya et al. 2017).

This work is supported by NASA Heliophysics Division through its Heliophysics Supporting Research (grant 80NSSC18K0066) and Guest Investigator (grant 80NSSC18K0068) programs, and by the Solar Terrestrial program of the National Science Foundation (grant AGS-1623844). Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. K.D. acknowledges helpful discussions with Michal Švanda.

A first step in computing flow kernels for the synthetic data is to obtain a model for the line widths and an empirical source function corresponding to these data. We begin from the azimuthally averaged power spectrum from the simulation (30 hr total). We fit the power spectrum of the synthetic data with a function of the form of Equation (53) from Birch et al. (2004), with the following assumptions: (1) the source correlation time is zero, this factor will be accounted for in the empirical source function, (2) instead of assuming a particular source depth we instead allow the source function (here, denoted s_n(k)) to be a free function in the fit, (3) we allow the mode frequencies to deviate from Model S (Christensen-Dalsgaard et al. 1996) frequencies, and (4) the damping rates γ_n(k) are also treated as free parameters. We carry out the fit in the range 0.1 rad Mm⁻¹ < k < 1.8 rad Mm⁻¹, using the range where ω/2π is between 2.5 mHz and 5.5 mHz and the horizontal phase speed is less then 60 km s⁻¹. To stabilize the fit, we parameterize the damping rates at each radial order and the real and imaginary parts of the source function at each radial order, and the deviation from the Model S mode frequencies at each radial order as sums of b-splines that are functions of horizontal wavenumber. We choose the number of b-splines for each radial order and physical quantity by hand, with the qualitative goal of a keeping the number of free parameters as small as possible while capturing the significant variations of the power spectrum of the synthetic data. In a few cases, the use of b-splines provides more freedom than is needed; in these cases, we used linear functions of k. Table 3 shows details of the choice of b-splines (or linear functions) for the fit. At each radial order, the knots of the b-splines are equally spaced between k_min and k_max for that radial order, with duplicate knots at the end points of the interval. The damping rates for radial orders of four or less are fit with functions of the form ${\gamma }_{n}(k)={{\rm{\Gamma }}}_{n,0}+{{\rm{\Gamma }}}_{n,1}{k}^{\alpha }$ with Γ_n,0, Γ_n,1, and α as free parameters and with a constant line width used for n = 5 and n = 6.

Though we allow the mode frequencies to vary in the fit, we use Model S stratification and eigenfunctions in the calculations of the Born-approximation kernels (Appendix B). Thus, the kernel calculations do not account for the difference between the mode frequencies in Model S and the mode frequencies in the synthetic data. This difference between the mode frequencies corresponds on changes in the horizontal phases speeds of the order of 2% and thus presumably an error in the kernels of the same order. Errors in the kernels of this amplitude imply errors in the forward-modeled travel times of less than 1 s, which is well below the error estimates for the travel times (Table 2). We thus expect that the errors in the inversions caused by the assumption of Model S stratification will be small compared to the errors caused by noise.

Figure 10 compares the model power spectrum (left) and the model power spectrum associated with the kernel calculation (right). For the part of the diagram with horizontal phase speed less then 55 km s⁻¹, there is reasonable qualitative agreement.

Zoom In Zoom Out Reset image size

Figure 11 shows two example slices through the power spectra at a constant horizontal wavenumber. In these slices, the result of the fit is shown along with the power spectra from the synthetic data and the model power spectrum from the kernel calculation. As discussed earlier, the kernel calculation is carried out using the Model S stratification, and so the resonance frequencies are slightly different than in the simulations.

Zoom In Zoom Out Reset image size

For the calculation of the kernels, we work in a coordinate system where ${\boldsymbol{r}}$ is horizontal position and z is height measured from the photosphere (z = 0). The background model is translation invariant and given by a plane-parallel version of Model S.

At each temporal frequency, the ingression H₋ at horizontal position r for a particular focus height z is related to the observed wavefield ϕ by

$$\begin{eqnarray}&&{H}_{-}({\boldsymbol{r}})=\int {G}_{-}^{P}({\boldsymbol{r}}^{\prime} -{\boldsymbol{r}})\phi ({\boldsymbol{r}}^{\prime} )\ d{\boldsymbol{r}}^{\prime} \end{eqnarray} \tag{ 7 }$$

where ${G}_{-}^{P}$ is the anti-causal Green's function multiplied by the appropriate pupil function P and data-analysis filter function. The Green's function, the pupil function, and the filter function all depend on the focus depth (see Section 3). For the sake of readability, we have also suppressed the notation showing that H₋, G₋, and ϕ are all functions of temporal frequency. The integral is taken over all horizontal positions where the pupil function is not zero. Equation (7) shows that the ingression is the result of filtering the wavefield with a non-axisymmetric filter (the 2D Fourier transform of ${G}_{-}^{P}$). The egression is related to the wavefield by an analogous equation, but with G₋ replaced by the causal Green's function G₊, and the pupil function replaced by its appropriate counterpart (e.g., for the NS travel-time difference, if P is the north pupil, then P' is the south pupil).

At each temporal frequency, The lateral-vantage ingression–egression covariance C at the horizontal position ${\boldsymbol{r}}=0$ is

$$\begin{eqnarray}&&C({\boldsymbol{r}})={H}_{-}^{* }({\boldsymbol{r}}){H}_{+}({\boldsymbol{r}}).\end{eqnarray} \tag{ 8 }$$

The ingression–egression covariance is a time–distance covariance at zero distance (recalling that the ingression and egression are filtered versions of the wavefield). We can use Equations (11)–(13), which already allow for arbitrary non-axisymmetric filters, from Birch & Gizon (2007) to compute the linear sensitivity of C to flows.

Travel-time shifts are measured from the ingression–egression covariance C as

$$\begin{eqnarray}&&\tau =\mathrm{Arg}[\bar{C}]/\bar{\omega },\end{eqnarray} \tag{ 9 }$$

where $\bar{C}$ is the sum of C over all frequency bins and $\bar{\omega }$ is the average frequency weighted by $| C|$. The perturbation δτ to the travel-time shift coming from a perturbation $\delta \bar{C}$ to the frequency-summed ingression–egression covariance is

$$\begin{eqnarray}&&\delta \tau =\mathrm{Im}[{\bar{C}}^{* }\delta \bar{C}]/[| \bar{C}{| }^{2}\bar{\omega }].\end{eqnarray} \tag{ 10 }$$

This equation shows that the linear sensitivity of the travel-time shift to flows is a linear combination of the kernels for the ingression–egression covariances (Equation (8). The resulting vector-valued kernels ${\boldsymbol{K}}$ satisfy

$$\begin{eqnarray}&&\delta \tau ({\boldsymbol{r}})=\int \int \int {\boldsymbol{K}}({\boldsymbol{r}}^{\prime} -{\boldsymbol{r}},z)\cdot {\boldsymbol{u}}({\boldsymbol{r}}^{\prime} ,z)\ d{\boldsymbol{r}}^{\prime} {dz}.\end{eqnarray} \tag{ 11 }$$

Figure 12 shows slices through the kernels ${K}_{\beta }^{\mathrm{we}}$ for the focus depth 3.97 Mm. The kernel ${K}_{x}^{\mathrm{we}}$ is symmetric in both x and y, the kernel ${K}_{y}^{\mathrm{we}}$ is anti-symmetric in both x and y, and the kernel ${K}_{z}^{\mathrm{we}}$ is anti-symmetric in x and symmetric in y. The ${K}_{x}^{\mathrm{we}}$ and ${K}_{z}^{\mathrm{we}}$ kernels have larger amplitudes than the ${K}_{y}^{\mathrm{we}}$ kernel; as expected, travel-time differences in the x-direction are mostly sensitive to flows in the x-direction and vertical flows. The depth dependence of the ${K}_{x}^{\mathrm{we}}$ is shown in more detail in Figure 13.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 14 shows slices through the kernels ${K}_{\beta }^{\mathrm{oi}}$ for the focus depth 3.97 Mm. The symmetries are different than for the case of travel-time difference in the x-direction. The kernel ${K}_{x}^{\mathrm{oi}}$ is anti-symmetric in x and symmetric in y. The kernel ${K}_{y}^{\mathrm{oi}}$ is symmetric in x and anti-symmetric in y. The kernel ${K}_{z}^{\mathrm{oi}}$ is cylindrically symmetric about the z-axis. The largest sensitivity is to vertical flows near the axis. The sign is such that upflows cause an increase in the OI time difference.

Zoom In Zoom Out Reset image size