© ESO, 2011

1. Introduction

The sub-surface of the Sun is optically thick, preventing us from directly observing the interior layers. Understanding the properties of the plasma in these regions has consequences for the theories of convection, stability of sunspots, the dynamics of stratified convection, and others. Most current knowledge about convection comes primarily from computational work (e.g., Vögler et al. 2005; Benson et al. 2006; Rempel et al. 2009). Helioseismic inversions of the sub-surface flows will play an important role in constraining these theories.

A powerful way of imaging the solar interior is via inferences gathered from studying the statistics of the acoustic and surface gravity waves at the surface. Solar pressure and surface gravity modes are generated randomly by the vigorous turbulence in the upper convection zone. These oscillations are observed in the solar photosphere by measuring Doppler shifts of photospheric absorption lines. Forward modelling allows us to relate anomalies (like flows, thermal hot/cold spots etc.) to changes in helioseismic observables.

The aim of helioseismic inversions is to reveal the structure of the subsurface flows (rotation, meridional circulation, convection), magnetic fields, and to measure deviations in the plasma state parameters (temperature, density, pressure) from a quiet Sun average. In this paper, we focus on travel times (Duvall et al. 1993), i.e., quantities that emerge from fits to cross correlations of observed signals. Time-distance helioseismology is used to measure and interpret changes in travel times of seismic waves caused by inhomogeneities in the structure of the Sun (see review by Gizon et al. 2010). In recent years, time-distance helioseismology has been used to invert for near-surface flows (e.g., Gizon et al. 2000; Duvall & Gizon 2000; Zhao & Kosovichev 2004; Jackiewicz et al. 2008), for flows beneath sunspots (e.g., Duvall et al. 1996; Zhao et al. 2001; Couvidat et al. 2006; Cameron et al. 2008; Gizon et al. 2009; Moradi et al. 2010) and flows in their vicinity (Gizon et al. 2000), study the rotational gradient at the base of the convection zone (e.g., Hanasoge & Duvall 2009), etc.

Helioseismic inversions are performed using two principal methods: the regularised least squares (RLS) and optimally localised averaging (OLA). The RLS method (in time-distance helioseismology used for the first time by Kosovichev 1996) seeks to find the models of the solar interior, which provide the best least-squares fit to the measured travel-time maps, while regularising the solution (e.g., by requiring the smooth solution). The OLA method was developed for geoseismology (Backus & Gilbert 1968, 1970). A form suitable for use in helioseismology was devised by Pijpers & Thompson (1992), who formulated the subtractive-OLA method. SOLA is based on explicitly constructed spatially confined averaging kernels by taking linear combination of sensitivity kernels, while simultaneously keeping the error magnification small. The resulting coefficients are then used to linearly combine the travel-time maps and obtain the estimate for structure and magnitude of solar plasma perturbations. A SOLA-type inversion is the principal method discussed in the current paper. The SOLA has been used in time-distance local helioseismology in the past by Jackiewicz et al. (2007, 2008) who demonstrated the ability of SOLA inversions to reveal the structure of the 3D internal flows. An efficient approach to solve fully consistent SOLA inversions was introduced by Jackiewicz et al. (2011). In this paper we focus on inversions for three-dimensional vector flows on supergranular scales in the near-surface layers of the solar convection zone.

1.1. Validating helioseismic inversions

Two approaches have been used so far to validate time-distance inversions. The first approach consists of generating synthetic travel-time maps by convolving a (known) arbitrary frozen flow field with travel-time sensitivity kernels and then testing the inversion method using these synthetic travel times (e.g., Kosovichev & Duvall 1997). The second approach (Zhao et al. 2007, 2010) is to use evolving realistic numerical simulations of three-dimensional radiative convection, where the helioseismic waves are naturally excited by the convection.

The first approach is convenient, but may not represent a realistic situation, in particular because the prescribed subsurface structures are often too idealised. The second approach is preferable, but is limited by computing resources: only simulations of the very near-surface regions of the Sun are available today. Both approaches have shown that inversions are generally able to retrieve the horizontal components of velocity at supergranulation scales in the quiet Sun. However, some problems have been reported. For example, vertical velocities have been measured with the opposite sign near the surface (Zhao et al. 2007) and the structure and sign of the flows around sunspots vary with the inversion method used (see, e.g., Gizon et al. 2009; Moradi et al. 2010).

The aim of this paper is to validate a particular implementation of SOLA inversions for time-distance helioseismology using a mixed approach: we take a snapshot from a large-box realistic simulation of solar convection (Ustyugov 2008) to generate realistic travel-time maps by convolution with Born travel-time sensitivity kernels. A realistic noise component is added to the travel times. The travel-time maps are then inverted using a multichannel SOLA inversion and compared with the known flows. This approach allows us to investigate various types of biases in the results of the inversion and to develop a robust procedure to minimise them. The most serious bias comes from the natural correlations among the components of flows induced by mass conservation that translate into a cross-talk between the components of the inverted flow (as we shall see, this is important when retrieving the vertical flow). Using realistic numerical simulations of convective flow velocities as input is very useful to set the acceptable level of random noise of the inverted flow velocities. The use of travel times with realistic noise properties is important in order to derive realistic estimates of the noise of the inverted flow velocities.

2. Synthetic travel times

To construct synthetic travel times, we use a realistic hydrodynamic simulation (Ustyugov 2008) of the solar convection. The computational domain is a box 20 Mm in depth and 60 Mm in each horizontal direction. The simulation provides us with a reasonably realistic description of flows in the upper convection zone v = (v_x,v_y,v_z). Throughout this paper the spatial coordinates are defined as ${x}\mathrm{=}\mathrm{\left(}{r}\mathit{,z}\mathrm{\right)}\mathit{,}$(1)where r is the horizontal position vector and z is the height.

We choose a snapshot (t = 500 min, see Fig. 16 in the on-line supplement) of the above mentioned simulation for our inversion tests. Following Duvall et al. (1997), we consider different types of travel-time measurements between a surface point at position r and a concentric annulus or quadrants in order to measure travel times sensitive to flows in inward-outward (denoted by “oi”), west-east (“we”), and north-south (“ns”) directions. Travel-time maps are denoted by τ^a(r), where the superscript a is an integer that uniquely refers to a particular combination of choices in the data analysis: a type of geometry (oi, we, or ns), annulus radius (from 7.3 Mm to 29.2 Mm every 1.46 Mm), and a wave filter (here, ridge filters for one of f, p₁, p₂, p₃, or p₄ modes). Thus, index a refers to one of M = 3 × 16 × 5 = 240 possibilities. Additional information describing the measurement procedure in time-distance helioseismology is given by Gizon & Birch (2005).

Travel-time maps (travel times as functions of position vector r) are generated by convolving the convection snapshot with sensitivity kernels according to ${\mathit{\tau }}^{\mathit{a}}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{\odot }}{K}\mathit{a}\mathrm{(}{r}\mathrm{\prime }\mathrm{-}{r}\mathit{,z}\mathrm{)}\mathrm{·}{v}\mathrm{\left(}{r}\mathrm{\prime }\mathit{,z}\mathrm{\right)}{\mathrm{d}}^{\mathrm{2}}{r}\mathrm{\prime }\hspace{0.17em}\mathrm{d}\mathit{z}\mathrm{+}{\mathit{n}}^{\mathit{a}}\mathrm{\left(}{r}\mathrm{\right)}\mathit{,}$(2)where ${K}{}^{\mathit{a}}\mathrm{=}\mathrm{\left(}{\mathit{K}}_{\mathit{x}}^{\mathit{a}}\mathit{,}{\mathit{K}}_{\mathit{y}}^{\mathit{a}}\mathit{,}{\mathit{K}}_{\mathit{z}}^{\mathit{a}}\mathrm{\right)}$ is a vector travel-time sensitivity kernel (see, e.g., Fig. 17 in the on-line supplement), v is the velocity vector of convecting flows, and the volume integral is taken over the Sun. The noise component is denoted by n^a. The original horizontal size of the simulation box of 60 Mm was too small for our purpose. Since the simulated velocities are periodic in the horizontal directions, we copied the simulation box 10 × 10 times.

Kernels are computed using the single-scattering Born approximation (Birch & Gizon 2007) and depend on eigenmodes of a background 1D standard solar model (Christensen-Dalsgaard et al. 1996). The kernels are invariant under horizontal translations; horizontal averages of these kernels for various oscillation modes are displayed in Fig. 18 in the on-line supplement to display their sensitivity in depth. All kernels used in this study have sensitivities only in upper-most 10 Mm of the convection zone.

Solar waves are excited by the action of turbulent convection in the Sun and therefore travel times are inherently noisy. Following Gizon & Birch (2004), we use a realistic noise covariance matrix (see Fig. 19 in the on-line supplement), ${\mathrm{\Lambda }}_{\mathit{ab}}\mathrm{(}{r}\mathrm{-}{r}\mathrm{\prime }\mathrm{)}\mathrm{=}\mathit{E}\mathrm{\left[}{\mathit{n}}^{\mathit{a}}\mathrm{\right(}{r}\mathrm{\left)}{\mathit{n}}^{\mathit{b}}\mathrm{\right(}{r}\mathrm{\prime }\mathrm{\left)}\mathrm{\right]}\mathit{,}$(3)to generate realizations of the noise.

Each travel-time map has 400 × 400 pixels with spatial sampling of 1.46 Mm, corresponding to the pixel size of a Michelson-Doppler-Imager full-disc image (MDI, Scherrer et al. 1995). In total, we generate 240 different travel-time maps (one for each index a, see, e.g., Fig. 1). These travel-time maps have spatial power spectra that are similar to observed travel-time power spectra based on MDI data, as demonstrated in Fig. 2. However, our synthetic travel times have weaker power at low wave-numbers due to the complete lack of signal from the simulated convective velocities below k   R_⊙ = 60 (only simulated noise is present).

We check to ensure that point-to-annulus sensitivity kernels and noise covariance matrices obey expected symmetries. Furthermore, the sensitivity kernels and covariance matrices are forced to decay smoothly to zero towards the edge of the computation box by multiplying with a smooth spatial function having zeros far from the region of interest. In addition, we make sure that the horizontal integral of each sensitivity kernel for v_z is zero at all depths. This is a consequence of symmetries associated with point-to-point kernels for v_z, which are insensitive to the mean value of vertical velocity (Birch & Gizon 2007).

3. Improved SOLA inversions

3.1. Real-space formulation

The inversion method we use here is an improved version of the subtractive optimally localised averages (SOLA) method described by Jackiewicz et al. (2007, 2008, 2011).

The SOLA algorithm describes how to optimally combine a given set of travel time measurements to infer the underlying properties of medium. In this paper, we wish to retrieve $\stackrel{}{{{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}}_{\mathrm{˜}}}\mathrm{\left(}{x}{}_{\mathrm{0}}\mathrm{\right)}$, i.e., an estimate of the α-component of the flow velocity in the neighbourhood of position x₀ = (r₀,z₀) in the solar interior. In practice, we search for a linear combination of the travel-time measurements, $\stackrel{}{{{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}}_{\mathrm{˜}}}\mathrm{\left(}{x}\mathrm{0}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i,a}}{\mathit{w}}_{\mathit{a}}^{\mathit{\alpha }}\mathrm{(}{r}\mathit{i}\mathrm{-}{r}\mathrm{0}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}{\mathit{\tau }}^{\mathit{a}}\mathrm{\left(}{r}\mathit{i}\mathrm{\right)}\mathit{,}$(4)where ${\mathit{w}}_{\mathit{a}}^{\mathit{\alpha }}\mathrm{\left(}{r}{}_{\mathit{i}}\mathrm{\right)}$ are inversion weights to be determined. Combining Eqs. (2) and (4), we have $\begin{array}{ccc}\stackrel{}{{{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}}_{\mathrm{˜}}}\mathrm{\left(}{x}\mathrm{0}\mathrm{\right)}& \mathrm{=}& \int \mathrm{\odot }\mathbf{𝒦}\mathit{\alpha }\mathrm{(}{r}\mathrm{-}{r}\mathrm{0}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}\mathrm{·}{v}\mathrm{\left(}{x}\mathrm{\right)}{\mathrm{d}}^{\mathrm{2}}{r}\hspace{0.17em}\mathrm{d}\mathit{z}\\ & & \end{array}$(5)where $\mathbf{𝒦}{}^{\mathit{\alpha }}\mathrm{=}\mathrm{\left(}{\mathrm{𝒦}}_{\mathit{x}}^{\mathit{\alpha }}\mathit{,}{\mathrm{𝒦}}_{\mathit{y}}^{\mathit{\alpha }}\mathit{,}{\mathrm{𝒦}}_{\mathit{z}}^{\mathit{\alpha }}\mathrm{\right)}$ is a vector averaging kernel with components ${\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{\equiv }\sum _{\mathit{i,a}}{\mathit{w}}_{\mathit{a}}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathit{i}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}{\mathit{K}}_{\mathit{\beta }}^{\mathit{a}}\mathrm{(}{r}\mathrm{-}{r}\mathit{i}\mathit{,z}\mathrm{)}\mathit{.}$(6)The expression d²r   dz means d³x.

The expectation value of inverted velocity $\stackrel{}{{{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}}_{\mathrm{˜}}}\mathrm{\left(}{x}{}_{\mathrm{0}}\mathrm{\right)}$ may be expanded as $\mathrm{E}\mathrm{\left[}\stackrel{}{{{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}}_{\mathrm{˜}}}\mathrm{\right]}\mathrm{=}{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}\mathrm{+}\sum _{\mathit{\beta }\ne \mathit{\alpha }}{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}{\mathrm{\left(}\mathit{\beta }\mathrm{\right)}}^{}\mathit{,}$(7)where ${\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}{\mathrm{\left(}\mathit{\beta }\mathrm{\right)}}^{}\mathrm{\left(}{x}\mathrm{0}\mathrm{\right)}\mathrm{\equiv }{\mathrm{\int }}_{\mathrm{\odot }}{\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{(}{r}\mathrm{-}{r}\mathrm{0}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}{\mathit{v}}_{\mathit{\beta }}\mathrm{\left(}{x}\mathrm{\right)}{\mathrm{d}}^{\mathrm{2}}{r}\hspace{0.17em}\mathrm{d}\mathit{z,}$(8)and ${\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}\mathrm{\equiv }{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}{\mathrm{\left(}\mathit{\alpha }\mathrm{\right)}}^{}$. The form of Eq. (7) will be useful later in this paper to optimize the inversion. Notice that the second term on the right side of (7) represents leakage of other flow components into the inverted one. This cross-talk may be a significant source of bias and needs to be studied.

We search for inversion weights w^α so that the vector averaging kernel resembles a user-supplied target . Because we wish to invert for the α-component of the flow, we choose the target with components: ${\mathrm{𝒯}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{(}{r}\mathrm{-}{r}\mathrm{0}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}\mathrm{=}\mathrm{𝒯}\mathrm{(}{r}\mathrm{-}{r}\mathrm{0}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}{\mathit{\delta }}_{\mathit{\alpha \beta }}\mathit{,}$(9)where is a function that peaks around x₀ = (r₀,z₀), diminishes rapidly away from that point, and has unit integral. δ_αβ is a Kronecker δ. Throughout this paper, we choose a simple Gaussian $\mathrm{𝒯}\mathrm{\left(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathrm{4}\ln \mathrm{2}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{{\mathit{\pi }}^{\mathrm{3}\mathit{/}\mathrm{2}}{\mathit{s}}_{\mathit{h}}^{\mathrm{2}}{\mathit{s}}_{\mathit{z}}}\exp \left[\mathrm{-}\frac{\mathrm{4}\ln \mathrm{2}}{{\mathit{s}}_{\mathit{h}}^{\mathrm{2}}}\mathrm{\parallel }{r}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{-}\frac{\mathrm{4}\ln \mathrm{2}}{{\mathit{s}}_{\mathit{z}}^{\mathrm{2}}}\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\right]\mathit{,}$(10)where s_h and s_z are full widths at half maximum (FWHM) of the Gaussian in horizontal and vertical directions respectively.

A successful inversion will return a value of ${\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}\mathrm{\left(}{x}{}_{\mathrm{0}}\mathrm{\right)}$ close to the target velocity ${\mathit{v}}_{\mathit{\alpha }}^{\mathrm{tgt}}\mathrm{\left(}{x}\mathrm{0}\mathrm{\right)}\mathrm{\equiv }{\mathrm{\int }}_{\mathrm{\odot }}\mathrm{𝒯}\mathrm{(}{r}\mathrm{-}{r}\mathrm{0}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}{\mathit{v}}_{\mathit{\alpha }}\mathrm{\left(}{x}\mathrm{\right)}{\mathrm{d}}^{\mathrm{2}}{r}\hspace{0.17em}\mathrm{d}\mathit{z,}$(11)and cross-talk velocities ${\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}{\mathrm{\left(}\mathit{\beta }\mathrm{\right)}}^{}\mathrm{\left(}{x}{}_{\mathrm{0}}\mathrm{\right)}$ close to zero for β ≠ α.

A list of the notations used for the various flow velocities referred to in the inversion procedure is given in Table 1.

The problem to be solved to obtain the w^α belongs to the class of constrained regularised optimisations. The terms to be regularised are the following:

• The misfit, i.e., how far the averaging kernel is from the desired target function, ${\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}{\mathit{w}}^{\mathit{\alpha }}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{\odot }}\mathrm{\left[}{\mathrm{𝒦}}_{\mathit{\alpha }}^{\mathit{\alpha }}\mathrm{\right(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}\mathrm{-}\mathrm{𝒯}\mathrm{(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}{\mathrm{d}}^{\mathrm{2}}{r}\hspace{0.17em}\mathrm{d}\mathit{z,}$(12)where ${\mathrm{𝒦}}_{\mathit{\alpha }}^{\mathit{\alpha }}$ is an implicit function of the weights w^α.

• The cross-talk quantifying the leakage of the signal from other flow components into the inverted one, $\mathit{XT}\mathrm{\left(}{\mathit{w}}^{\mathit{\alpha }}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\beta }\ne \mathit{\alpha }}{\mathrm{\int }}_{\mathrm{\odot }}\mathrm{\left[}{\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{\right(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}{\mathrm{d}}^{\mathrm{2}}{r}\hspace{0.17em}\mathrm{d}\mathit{z}\mathit{.}$(13)

• The variance of $\stackrel{}{{{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}}_{\mathrm{˜}}}$, which corresponds to the root-mean-square of noise-related fluctuations in inferred velocities. We refer to this quantity as predicted error hereafter, ${\mathit{\sigma }}_{\mathit{\alpha }}^{\mathrm{2}}\mathrm{\left(}{\mathit{w}}^{\mathit{\alpha }}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i,j,a,b}}{\mathit{w}}_{\mathit{a}}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathit{i}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}{\mathrm{\Lambda }}_{\mathit{ab}}\mathrm{(}{r}\mathit{i}\mathrm{-}{r}\mathit{j}\mathrm{)}{\mathit{w}}_{\mathit{b}}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathit{j}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathit{.}$(14)

• The ad hoc term quantifying the spread of the weights in space around r₀, $\mathit{S}\mathrm{\left(}{\mathit{w}}^{\mathit{\alpha }}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{a,i}}\mathrm{\left[}{\mathit{w}}_{\mathit{a}}^{\mathit{\alpha }}\mathrm{\right(}{r}\mathit{i}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}\mathit{.}$(15)The regularisation based on S ensures that weights decrease towards the edge of the horizontal domain in order to enforce the spatial locality of the inversion and to prevent the weights from oscillating in the spatial domain.

In practice, we search for the weights w^α that minimise the cost function $\begin{array}{ccc}{\mathit{F}}^{\mathit{\alpha }}\mathrm{\left(}{\mathit{w}}^{\mathit{\alpha }}\mathrm{\right)}& \mathrm{=}& \int \mathrm{\odot }\mathrm{\left[}{\mathrm{𝒦}}_{\mathit{\alpha }}^{\mathit{\alpha }}\mathrm{\right(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}\mathrm{-}\mathrm{𝒯}\mathrm{(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}{\mathrm{d}}^{\mathrm{2}}{r}\hspace{0.17em}\mathrm{d}\mathit{z}\\ & & \mathrm{+}\hspace{0.17em}\mathit{\nu }\hspace{0.17em}\sum _{\mathit{\beta }\ne \mathit{\alpha }}{\mathrm{\int }}_{\mathrm{\odot }}\mathrm{\left[}{\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{\right(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}{\mathrm{d}}^{\mathrm{2}}{r}\hspace{0.17em}\mathrm{d}\mathit{z}\mathrm{+}\mathit{ϵ}\sum _{\mathit{a,i}}\mathrm{[}{\mathit{w}}_{\mathit{a}}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathit{i}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}{\mathrm{]}}^{\mathrm{2}}\\ & & \mathrm{+}\hspace{0.17em}\mathit{\mu }\sum _{\mathit{i,j,a,b}}{\mathit{w}}_{\mathit{a}}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathit{i}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}{\mathrm{\Lambda }}_{\mathit{ab}}\mathrm{(}{r}\mathit{i}\mathrm{-}{r}\mathit{j}\mathrm{)}{\mathit{w}}_{\mathit{b}}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathit{j}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathit{.}\end{array}$(16)Trade-off parameters μ, ν, and ϵ balance these terms. The strategy to be employed in order to set these values will be described in Sect. 3.3.

The inversion is subject to constraints ${\mathrm{\int }}_{\mathrm{\odot }}{\mathrm{d}}^{\mathrm{2}}{r}\hspace{0.17em}\mathrm{d}\mathit{z}\hspace{0.17em}{\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}{\mathit{\delta }}_{\mathit{\alpha \beta }}\mathrm{forall}\mathit{\beta ,}$(17)in order to scale the amplitude of the inverted flow $\begin{array}{}{}_{\mathbf{˜}}\\ {v}\end{array}{}^{\mathrm{inv}}$ appropriately.

By taking derivative of (16) with constraint (17) added with respect to the weights, the problem can be cast into a linear inverse problem, as explained by Pijpers & Thompson (1992) or Jackiewicz et al. (2011). In real space, the matrix to be inverted has (N²M + P)² ≃ (10⁷)² elements, where N = 200 is the number of grid points in one horizontal direction, P = 3 is the number of physical unknowns (three velocity components), and M = 240 as already defined earlier.

3.2. Fourier-space (multichannel) formulation

The full problem written in real space is intractable to be solved using nowadays computers, because the matrix to be inverted is too large. Jackiewicz et al. (2011) found a solution to this problem by transforming to spatial Fourier space, where the inverse problem decouples as a consequence of the horizontal translation invariance of the sensitivity kernels. Thus, instead of a big linear inverse problem (inverting matrix having (10⁷)² elements), we solve 40   000 small linear inverse problems (inverting matrices having 240² elements) in wave-vector space. This approach is called a multichannel inversion (Jensen et al. 1998). Here we summarise results of Jackiewicz et al. (2011) for the sake of completeness.

Following Jackiewicz et al. (2011), we use the following definition of the Fourier transform, such that any function f(r) and its 2D spatial Fourier transform are related according to $\begin{array}{ccc}\mathit{f}\mathrm{\left(}{r}\mathrm{\right)}& \mathrm{=}& {\mathit{h}}_{\mathit{k}}^{\mathrm{2}}\sum _{{k}}\mathit{f̅}\mathrm{\left(}{k}\mathrm{\right)}\exp \left(\mathrm{i}{k}\mathrm{·}{r}\right)\mathit{,}\\ \mathit{f̅}\mathrm{\left(}{k}\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{h}}_{\mathit{x}}^{\mathrm{2}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}}\sum _{{r}}\mathit{f}\mathrm{\left(}{r}\mathrm{\right)}\exp \left(\mathrm{-}\mathrm{i}{k}\mathrm{·}{r}\right)\mathit{,}\end{array}$where h_x = 1.46 Mm and h_k = 0.022 rad/Mm are the grid spacings in the real and Fourier domains respectively.

For each non-vanishing wave-vector k = (k_x,k_y), the vector of weights $\mathit{W}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}\mathrm{\left[}\mathit{w̅}\begin{array}{}\mathit{\alpha }\\ \mathrm{1}\end{array}\mathrm{\right(}{k}\mathrm{\left)}\mathit{w̅}\begin{array}{}\mathit{\alpha }\\ \mathrm{2}\end{array}\mathrm{\right(}{k}\mathrm{\left)}\mathit{...}\mathit{w̅}\begin{array}{}\mathit{\alpha }\\ \mathit{M}\end{array}\mathrm{\right(}{k}\mathrm{)}{\mathrm{]}}^{\mathrm{T}}$ is the solution to the matrix equation, ${\mathit{h}}_{\mathit{k}}^{\mathrm{2}}{\mathit{N}}^{\mathrm{2}}\hspace{0.17em}\mathit{A}\mathrm{\left(}{k}\mathrm{\right)}\mathit{W}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}\mathit{T}\mathrm{\left(}{k}\mathrm{\right)} \mathrm{for}{k}\ne \mathbf{0}\mathit{.}$(20)Each matrix A(k) has M × M elements A_ab given by $\begin{array}{ccc}{\mathit{A}}_{\mathit{ab}}\mathrm{\left(}{k}\mathrm{\right)}& \mathrm{=}& \mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}\underset{\mathrm{-}\mathrm{\infty }}{\overset{\mathrm{+}\mathrm{\infty }}{\mathrm{\int }}}\mathrm{d}\mathit{z}\mathit{K̅}\begin{array}{}\mathit{a}\mathrm{\ast }\\ \mathit{\alpha }\end{array}\mathrm{\right(}{k}\mathit{,z}\mathrm{\left)}\mathit{K̅}\begin{array}{}\mathit{b}\\ \mathit{\alpha }\end{array}\mathrm{\right(}{k}\mathit{,z}\mathrm{)}\\ & & \mathrm{+}\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}\mathit{\nu }\sum _{\mathit{\beta }\ne \mathit{\alpha }}\underset{\mathrm{-}\mathrm{\infty }}{\overset{\mathrm{+}\mathrm{\infty }}{\mathrm{\int }}}\mathrm{d}\mathit{z}\mathit{K̅}\begin{array}{}\mathit{a}\mathrm{\ast }\\ \mathit{\beta }\end{array}\mathrm{\right(}{k}\mathit{,z}\mathrm{\left)}\mathit{K̅}\begin{array}{}\mathit{b}\\ \mathit{\beta }\end{array}\mathrm{\right(}{k}\mathit{,z}\mathrm{)}\\ & & \end{array}$(21)The vector has M elements given by $\mathit{t̅}\mathit{a}\mathrm{\left(}{k}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}\underset{\mathrm{-}\mathrm{\infty }}{\overset{\mathrm{+}\mathrm{\infty }}{\mathrm{\int }}}\mathrm{d}\mathit{z}\hspace{0.17em}\mathit{K̅}\begin{array}{}\mathit{a}\mathrm{\ast }\\ \mathit{\alpha }\end{array}\mathrm{(}{k}\mathit{,z}\mathrm{\left)}\mathrm{𝒯̅}\mathrm{\right(}{k}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{)}\mathit{.}$(22)In addition, the matrix equation for the case k = 0 is $\left[\begin{array}{c}\\ {\mathit{h}}_{\mathit{k}}^{\mathrm{4}}{\mathit{N}}^{\mathrm{2}}\mathit{A}\mathrm{\left(}\mathbf{0}\mathrm{\right)}& \mathit{C}\\ {\mathit{C}}^{\mathrm{T}}& \mathrm{0}\end{array}\right]\left[\begin{array}{c}\\ \mathit{W}\mathrm{\left(}\mathbf{0}\mathrm{\right)}\\ \mathit{L}\end{array}\right]\mathrm{=}\left[\begin{array}{c}\\ {\mathit{h}}_{\mathit{k}}^{\mathrm{2}}\hspace{0.17em}\mathit{T}\mathrm{\left(}\mathbf{0}\mathrm{\right)}\\ {\mathit{U}}^{\mathit{\alpha }}\mathit{/}\mathrm{\left(}{\mathit{h}}_{\mathit{k}}^{\mathrm{2}}{\mathit{N}}^{\mathrm{2}}\mathrm{\right)}\end{array}\right]\mathit{,}$(23)where C is an M × P matrix with elements ${\mathit{C}}_{\mathit{a\beta }}\mathrm{=}{\mathrm{\int }}_{\mathrm{\odot }}{\mathrm{d}}^{\mathrm{2}}{r}\hspace{0.17em}\mathrm{d}\mathit{z}\hspace{0.17em}{\mathit{K}}_{\mathit{\beta }}^{\mathit{a}}\mathrm{\left(}{r}\mathit{,z}\mathrm{\right)}\mathit{,}$(24)L is a 1 × P vector of Lagrange multipliers, and the vector U^α is a 1 × P unit vector with components ${\mathit{u}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{=}{\mathit{\delta }}_{\mathit{\alpha \beta }}$.

Solutions to the N − 1 Eq. (20) and the Eq. (23) give Fourier components of the weights, $\mathit{w̅}\begin{array}{}\mathit{\alpha }\\ \mathit{a}\end{array}\mathrm{\left(}{k}\mathrm{\right)}$. By taking an inverse Fourier transform we obtain ${\mathit{w}}_{\mathit{a}}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathrm{\right)}$. Equation (4) then gives an estimate of the α-component of the flow, $\stackrel{}{{{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}}_{\mathrm{˜}}}$.

The difference between the above equations and those in Jackiewicz et al. (2011) arises from additional constraint terms (13) and (15).

3.3. Picking trade-off parameters

Trade-off parameters μ, ν, and ϵ control the balance between various terms in the cost-function (16).

In practice, during the inversion for each flow component, we compute a grid of solutions by varying all three trade-off parameters. For each solution on this grid, we compute the misfit (12), predicted error of results (14), amount of cross-talk (13), and the spatial power of weights (15). All four quantities can be computed in spatial Fourier domain.

Standard optimisation methods to achieve optimality, such as L-curve analysis (see, e.g., Hansen 1998; Jackiewicz et al. 2008), are not particularly useful to our problem. The elbow of the L-curve, which is considered an optimal point in parameter space, was located where the predicted noise level was very large. We therefore developed our own strategy of selecting values for the trade-off parameters.

We start with parameter μ, which controls the trade-off between the misfit and random error. We select μ such that the noise level of the inverted quantity is less than some target value, chosen by the user, depending on the problem at hand. For example, one may choose a target noise level of 20 m/s to invert for horizontal flows in supergranules. A much lower noise level will be required to invert for the vertical component of the flow (Duvall & Birch 2010). The example noise levels σ_α are given in Table 2, where we give the predicted inversion error for $\stackrel{}{{{\mathit{v}}_{\mathit{x}}^{\mathrm{inv}}}_{\mathrm{˜}}}$ and $\stackrel{}{{{\mathit{v}}_{\mathit{z}}^{\mathrm{inv}}}_{\mathrm{˜}}}$ and the root-mean-square of the flow averaged with the target function (${\mathit{v}}_{\mathit{x}}^{\mathrm{tgt}}$ and ${\mathit{v}}_{\mathit{z}}^{\mathrm{tgt}}$) and averaged with the resulting averaging kernel (${\mathit{v}}_{\mathit{x}}^{\mathrm{inv}}$ and ${\mathit{v}}_{\mathit{z}}^{\mathrm{inv}}$). It may be that this selection of μ leads to an averaging kernel which does not resemble the desired target function at all, in which case one would have to allow for a coarser spatial resolution.

For the given μ, we then choose ν so that the degree of cross-talk is less than 10^-5 Mm^-3. We observe that this constraint places the following upper bound on the magnitude of the cross-talk terms: $\frac{\max \left|\sum _{\mathit{\beta }\ne \mathit{\alpha }}{\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\right|}{\max \left|{\mathrm{𝒦}}_{\mathit{\alpha }}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\right|}\mathit{<}\mathrm{0.05.}$(25)The parameter ϵ controls the degree of spatial confinement of the inversion weights and has an impact on the misfit and error as well. If ϵ is too large, the weights will be highly localised in space around the central point, but a large misfit will result.

For the particular cases that we study here, we find that localisation of weights is accomplished when S < 2 × 10^-3   km²   s^-2. The impact of this regularisation term is illustrated in Fig. 3. For the case shown, weights are confined to within a disc of radius  ~50 Mm.

3.4. The code

The code implementing the above described procedure is written in Matlab. Matlab provides a compromise between computational efficiency and the availability of higher-level software constructs that make the code lucid, modular, and easy to modify. Matlab-based code may be compiled into binary-executable form suitable for pipeline pseudo-automatic processing. The input is a text file, in which the user specifies all relevant parameters including the kernels to be used in the inversion and the set of trade-off parameters to be investigated. This allows for the development of user interfaces such as the web-based or graphical interface, which can serve as the user front-end of the inversion code.

The code is parallel and scales linearly with number of parallel jobs. The execution is fast, e.g., a set of inversions for one flow component involving 240 kernels, 200  ×  200 spatial points, and a grid of 200 trade-off parameter values takes around 6 h using 48 Opteron-2.3 GHz CPUs. When many CPUs are involved, the extensive input-output load becomes a bottle-neck and affects the total processing time.

4. Inversion for horizontal flow

We applied the above-described method to known synthetic travel times in order to validate the set-up and the performance of the inversion procedure. We focus first on the horizontal flow components, i.e., α = x or α = y. We compare two principal quantities: (1) the flow map that is obtained by convolving the known velocity field with the target function (${\mathit{v}}_{\mathit{\alpha }}^{\mathrm{tgt}}$), i.e., the best-case inversion scenario, and (2) the flow map actually resulting from the inversion ($\stackrel{}{{{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}}_{\mathrm{˜}}}$). In this manner, we may investigate in detail different sources of bias in the results.

The simulation convolved with the target function gives us ${\mathit{v}}_{\mathit{\alpha }}^{\mathrm{tgt}}$, from which we estimate the expected magnitude of the flow that we want to invert for (see Table 2). This places limits on the required noise level of the inversion so that the results have signal-to-noise ratios larger than 1. Requirements on the targeted error level fix the trade-off parameter μ.

The choice of the target function depends on the discretion and needs of the user. Here we focus on layers in the top few Mm, in particular on depths  − z₀ = 1, 3.5, and 5.5 Mm. This set-up was selected because similar flow inversions were also performed and discussed by Jackiewicz et al. (2008), which makes it possible to compare the results of both methods. As discussed in the preceding sections, the outcome of the inversion is a set of weights ${\mathit{w}}_{\mathit{a}}^{\mathit{\alpha }}\mathrm{\left(}{r}\mathrm{\right)}$, which are used to combine the travel-time maps in order to obtain estimates of velocity $\stackrel{}{{{\mathit{v}}_{\mathit{x}}^{\mathrm{inv}}}_{\mathrm{˜}}}\mathrm{\left(}{r}\mathrm{\right)}$ and $\stackrel{}{{{\mathit{v}}_{\mathit{y}}^{\mathrm{inv}}}_{\mathrm{˜}}}\mathrm{\left(}{r}\mathrm{\right)}$. Example weights for $\stackrel{}{{{\mathit{v}}_{\mathit{x}}^{\mathrm{inv}}}_{\mathrm{˜}}}$ are displayed in Fig. 3.

All components of the averaging kernels ${\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{x}}\mathrm{\left(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}$ are shown in Figs. 20–22. We also show comparisons between inversions when the cross-talk is minimised and not minimised. The minimisation of the cross-talk with ν = 100 is very efficient. However, for depths of 3.5 and 5.5 Mm the minimisation of the cross-talk introduces some small artefacts in the inversion averaging kernel ${\mathrm{𝒦}}_{\mathit{x}}^{\mathit{x}}$, which are a small price to pay.

Figure 4 shows vertical cuts through the target functions and the averaging kernels when the cross-talk is minimised. Except for the target depth of 5.5 Mm, the depth dependences of the averaging kernels ${\mathrm{𝒦}}_{\mathit{x}}^{\mathit{x}}$ resemble that of the target functions. The near-surface inversion at 1 Mm is dominated by the f-mode. The inversion at 3.5 Mm has equal contributions from f and p₁, with opposite signs, as shown in Fig. 23 in the on-line supplement, where we plot the contribution of individual modes to the averaging kernels ${\mathrm{𝒦}}_{\mathit{x}}^{\mathit{x}}$ as a function of depth.

The validation of the v_x inversion is demonstrated in Fig. 5. Here we plot the desired ${\mathit{v}}_{\mathit{x}}^{\mathrm{tgt}}$ at three different depths and the inverted ${\mathit{v}}_{\mathit{x}}^{\mathrm{inv}}$ without noise contributions. These two are very close for the depths 1 and 3.5 Mm. The differences between ${\mathit{v}}_{\mathit{x}}^{\mathrm{tgt}}$ and ${\mathit{v}}_{\mathit{x}}^{\mathrm{inv}}$ at depth 5.5 Mm are caused by an imperfect averaging kernel. When the random noise is added to the solution (bottom row of Fig. 5), we see that the inversions for v_x at 1 Mm and 3.5 Mm are still very good, while the inverted v_x at depth 5.5 Mm is dominated by noise. Minimising cross-talk does not improve the quality of the solution in this case.

In order to quantify the inversion biases, we compare directly the expected and inverted values for a set of spatial locations. We select points separated by 7.5 Mm in each horizontal coordinate. This sampling interval is equal to the half of the horizontal FWHM of the target function (thus the points are somewhat independent). In addition, to avoid possible edge effects, we cut the outer part of the horizontal plane so that only the central 200  ×  200 pixel patch is kept. The scatter plots comparing various inversion components are displayed in Fig. 6 for the depth of 1 Mm, the plots for the depth of 3.5 Mm are qualitatively similar. We estimate that for the depth of 1 Mm, the imperfect averaging kernel leads to an average underestimation of the horizontal flow components by some 20% (Fig. 6 left). Furthermore we note that in the case of the inversion for the horizontal flow components, the bias caused by the cross-talk is not important and the random noise level corresponds to the predicted value (Fig. 6 right). The results for the depth 3.5 Mm are similar, while the results for the depth 5.5 Mm are dominated by the random noise.

The noise domination of v_x inversion at depth of 5.5 Mm is evident from the plot of the azimuthally averaged power spectra of the random noise and of the signal as a function of k   R_⊙ (Fig. 7 right). For comparison, similar plots for the depths 1 Mm and 3.5 Mm are given also in Fig. 7, where the power spectrum of the signal is well above the random noise near supergranular spatial scales. The decrease in power of the signal at low kR_⊙ in Fig. 7 is because the convection simulation does not contain these scales.

We showed that it is possible to retrieve v_x and v_y in the top 3.5 Mm without noticeable bias and within the predicted noise level of  ~25 m/s (for observing time T = 4 days) with 240 different travel-time measurements for ridges f to p₄. In Table 3, we summarise our findings. We state the statistical quantities (correlation coefficient and the slope of the linear fit) comparing $\stackrel{}{{{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}}_{\mathrm{˜}}}$ and the corresponding ${\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}$. We also show the signal-to-noise ratio of the inverted $\stackrel{}{{{\mathit{v}}_{\mathit{\alpha }}^{\mathrm{inv}}}_{\mathrm{˜}}}$.

The cross-talk is unimportant in inversions for the horizontal flow components. The cross-talk could come only from v_z, which is weak. We estimate that the influence of the cross-talk is less than 5%. The inversion at 5.5 Mm is already dominated by noise. If we drop the requirement on horizontal resolution, inversions at a depth of 5.5 Mm are also possible (e.g., with a FWHM of s_h = 25 Mm).

It is possible to perform inversions with less temporal averaging, but to obtain a reasonable signal-to-noise ratio, we would have to relax the demand on the match between the averaging kernel and the target function. Typically, this leads to side-lobes in the averaging kernel (especially in the z-direction), which may make interpretation of the results more difficult. Another possibility would be to include a larger number of independent travel-time measurements (and therefore more sensitivity kernels). These issues are being worked on.

5. Inversion for vertical flow

5.1. Specificity of inversions for vertical flow

Inversions for v_z require a different methodology because sensitivity kernels ${\mathit{K}}_{\mathit{z}}^{\mathit{a}}$ have zero horizontal integrals at each depth (and therefore zero total integral) implying (see Eq.(6)) $\underset{\mathrm{-}\mathrm{\infty }}{\overset{\mathrm{+}\mathrm{\infty }}{\mathrm{\int }}}{\mathrm{d}}^{\mathrm{2}}{r}\hspace{0.17em}{\mathrm{𝒦}}_{\mathit{z}}^{\mathit{z}}\mathrm{\left(}{r}\mathit{,}{\mathit{z}}^{\mathrm{\prime }}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\mathrm{forall}{\mathit{z}}^{\mathrm{\prime }}\mathit{.}$(26)Consequently, with measurements discussed here, it is impossible to retrieve horizontal average  ⟨ v_z ⟩  of vertical flow and we may only invert for fluctuations v_z −  ⟨ v_z ⟩ . Equation (26) implies

${\mathrm{\int }}_{\mathrm{\odot }}{\mathrm{d}}^{\mathrm{2}}{r}\hspace{0.17em}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\hspace{0.17em}{\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{z}}\mathrm{\left(}{r}\mathit{,}{\mathit{z}}^{\mathrm{\prime }}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathrm{for}\mathit{\beta }\mathrm{=}\mathit{z}$(27)and thus Eq. (17) cannot be written for α = z. This implies in turn that the matrix Eq. (23) for k = 0 may not be written either.

5.2. Ignoring k=0

A solution to this issue is simply to replace Eq. (23) by $\mathit{w̅}\begin{array}{}\mathit{z}\\ \mathit{a}\end{array}\mathrm{\left(}\mathbf{0}\mathrm{\right)}\mathrm{\equiv }\mathrm{0}\mathrm{forall}\mathit{a}\mathit{.}$(28)Since the averaging kernel ${\mathrm{𝒦}}_{\mathit{z}}^{\mathit{z}}$ must have zero integral, the averaging kernel will be offset by a small negative constant away from the central peak. Luckily, regularization by the term (15) will ensure that the averaging kernel will drop to zero at the edge of the inversion box. However, this solution is not quite satisfactory, because we do not have very much control over this extended negative surrounding sidelobe.

5.3. Target function with zero mean

It is more elegant to select a target function with vanishing horizontal integral at each depth. If this target function is to peak around x₀ = (r₀,z₀), then this peak must be compensated by negative side-lobes in horizontal directions. These side-lobes have to be constructed in such a way that they do not lead to significant biases. We suggest to choose $\mathrm{𝒯}\mathrm{\left(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathit{H}\mathrm{\left(}{r}\mathrm{\right)}\frac{\sqrt{\mathrm{4}\ln \mathrm{2}}}{\sqrt{\mathit{\pi }}{\mathit{s}}_{\mathit{z}}}\exp \left[\mathrm{-}\frac{\mathrm{4}\ln \mathrm{2}}{{\mathit{s}}_{\mathit{z}}^{\mathrm{2}}}\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}\right]\mathit{,}$(29)where $\begin{array}{ccc}\mathit{H}\mathrm{\left(}{r}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{4}\ln \mathrm{2}}{\mathit{\pi }{\mathit{s}}_{\mathit{h}}^{\mathrm{2}}}\exp \left[\mathrm{-}\frac{\mathrm{4}\ln \mathrm{2}}{{\mathit{s}}_{\mathit{h}}^{\mathrm{2}}}\mathrm{\parallel }{r}{\mathrm{\parallel }}^{\mathrm{2}}\right]\\ & & \end{array}$(30)The constant $\mathit{c}\mathrm{=}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}{\mathit{n}}^{\mathrm{2}}}\right)\mathrm{+}\sqrt{\frac{\mathit{\pi }}{\mathrm{2}{\mathit{n}}^{\mathrm{2}}}}\left[\mathrm{1}\mathrm{+}\mathrm{erf}\left(\frac{\mathrm{1}}{\mathit{n}\sqrt{\mathrm{2}}}\right)\right]$(31)ensures that the horizontal integral of is zero. The horizontal part of the target function is constructed from a Gaussian peaked at r₀ = 0, from which a wide surrounding Gaussian annulus is subtracted. The free parameter n balances the width of a side-lobe compared to the width of the central peak.

The inversion is performed for each k ≠ 0 with this new target function, together with $\mathit{w̅}\begin{array}{}\mathit{z}\\ \mathit{a}\end{array}\mathrm{(}{k}\mathrm{=}\mathbf{0}\mathrm{)}\mathrm{=}\mathrm{0}$.

5.4. Both target functions provide similar answers

We compared the results obtained using the two target functions proposed above and found two conclusions.

The two solutions described above provide results which are very close (see, e.g., Fig. 24). However, the solution with the target function having zero horizontal average is elegant and provides more control over the solution averaging kernel.

By comparing the results obtained with various sizes of the negative annulus surrounding the central peak we found that the bias in the inverted flow caused by this negative sidelobe is negligible if n > 3 in Eq. (30).

The results presented in the following sections are obtained using the first formalism.

5.5. Validation of vertical flow inversion

As is evident already, inversions for vertical flow are not as seamless as in the case of horizontal components. This is mostly because the vertical flow is much weaker on supergranular scales than the horizontal flow. As shown in Table 2, the expected rms of the vertical flow in the top layers is of the order of 5 m s^-1. Therefore, the predicted noise of the inverted v_z has to be set to a much smaller value than for the horizontal velocities. As a result, the match between the desired target function and the resulting averaging kernel will have to be poorer.

Furthermore, it is absolutely crucial to minimise the cross-talk to avoid the leakage of the large horizontal velocities into the small inferred vertical velocity. Even an apparently negligible cross-talk averaging kernel could in the end cause a significant bias in the results. Minimising the cross-talk is especially important because of the natural correlations between the vertical and horizontal flow components in the mass-conserving flow of the supergranules (see Fig. 8). In the upper layers, horizontal outflows are associated with upflows. The structure of the cross-talk averaging kernel ${\mathrm{𝒦}}_{\mathit{x}}^{\mathit{z}}$ as shown in the top half of Fig. 10 in the case when cross-talk is not minimised then implies a negative bias in v_z due to ${\mathit{v}}_{\mathit{z}}^{\mathrm{inv}}{\mathrm{\left(}\mathit{x}\mathrm{\right)}}^{}$ and ${\mathit{v}}_{\mathit{z}}^{\mathrm{inv}}{\mathrm{\left(}\mathit{y}\mathrm{\right)}}^{}$.

Figure 9 shows vertical cuts through the target functions and the averaging kernels when the cross-talk is minimised. The averaging kernels ${\mathrm{𝒦}}_{\mathit{z}}^{\mathit{z}}$ for the inversions at depth 3.5 and 5.5 Mm have sidelobes towards the surface. All components of the averaging kernel ${\mathrm{𝒦}}_{\mathit{\beta }}^{\mathit{z}}\mathrm{\left(}{r}\mathit{,z}\mathrm{;}{\mathit{z}}_{\mathrm{0}}\mathrm{\right)}$ for v_z are displayed in Fig. 10 (and Figs. 25, 26), again comparing the cases when the cross-talk is and is not minimised. The action of the cross-talk minimisation term is very efficient.

The validation of the v_z inversion is demonstrated in Fig. 11. Here we plot ${\mathit{v}}_{\mathit{z}}^{\mathrm{tgt}}$ at three different depths and the inverted ${\mathit{v}}_{\mathit{z}}^{\mathrm{inv}}$ without noise contributions. These two are very similiar at the depth of 1 Mm. Small differences between ${\mathit{v}}_{\mathit{z}}^{\mathrm{tgt}}$ and ${\mathit{v}}_{\mathit{z}}^{\mathrm{inv}}$ caused by an imperfect averaging kernel are visible at the depth of 3.5 Mm, become more significant at depth 5.5 Mm. When random noise is added to the solution (bottom row of Fig. 11), we see that the inversion for vertical flow is possible at 1 Mm depth only if cross-talk is minimised. Inversions at depths 3.5 and 5.5 Mm are buried in random noise, where minimising cross-talk does not help. The magnitude of cross-talk at a depth 1 Mm is much larger than that of the vertical flow; further, cross-talk is highly anti-correlated with the vertical flow, which makes the correlation coefficient between ${\mathit{v}}_{\mathit{z}}^{\mathrm{tgt}}$ and $\stackrel{}{{{\mathit{v}}_{\mathit{z}}^{\mathrm{inv}}}_{\mathrm{˜}}}$ close to  − 1. Note that a similarly high anti-correlation was measured by Zhao et al. (2007).

All components ${\mathit{v}}_{\mathit{z}}^{\mathrm{inv}}{\mathrm{\left(}\mathit{\beta }\mathrm{\right)}}^{}$ and the random noise component in the inversion for v_z at 1 Mm depth are displayed in Fig. 12. We see in the top row that the leakage of the horizontal components (${\mathit{v}}_{\mathit{z}}^{\mathrm{inv}}{\mathrm{\left(}\mathit{x}\mathrm{\right)}}^{}$ and ${\mathit{v}}_{\mathit{z}}^{\mathrm{inv}}{\mathrm{\left(}\mathit{y}\mathrm{\right)}}^{}$) covers up completely the weak signal of the vertical flow when the cross-talk is not minimised.

Biases in the v_z inversion at a depth of 1 Mm may be quantified by directly comparing expected and inverted values for a set of spatial locations in the horizontal plane (Fig. 13), as in the case of the v_x inversion (Sect. 4). The imperfect averaging kernel in this case does not cause any significant bias (Fig. 13 left). The v_z inversion is not sensitive to the horizontally averaged vertical flow  ⟨ v_z ⟩ , which is therefore subtracted from ${\mathit{v}}_{\mathit{z}}^{\mathrm{tgt}}$ in the corresponding plot. The effect of minimising cross-talk is shown in Fig. 13 left. The inverted $\stackrel{}{{{\mathit{v}}_{\mathit{z}}^{\mathrm{inv}}}_{\mathrm{˜}}}$ is anti-correlated with ${\mathit{v}}_{\mathit{z}}^{\mathrm{inv}}$ expected from the inversion due to the leakage of the horizontal mass-conserving flow components into the vertical one.

Meaningful inversions for the vertical flow are possible in shallow near-subsurface layers of the convection zone only when the cross-talk between the vertical and the horizontal components is minimised. Vertical flow inversions on supergranular horizontal scales at depths greater than  ~1 Mm performed using f to p₄ modes and travel time maps averaged over 4 days are dominated by random noise and the signal of the vertical flow may not therefore be inferred at all (see summary in Table 3). This fact is also demonstrated in Fig. 14, where we plot power spectra of individual inversion components as a function of the spatial scale. At depths larger than 1 Mm, the power of the signal is much less than the power of the random noise regardless of the spatial scale.

6. Beating the noise: statistical averaging

As demonstrated in previous sections, it is very difficult to have a meaningful inversion for solar flows on supergranular scales averaged over a short time, even in the shallow near-subsurface layers of the solar convection zone. The signal is overwhelmed by random noise. It is not feasible to average over longer times, because the expected lifetime of convection on supergranular scales is on the order of a day. By averaging over longer time, a significant portion of the scientifically useful information would be lost.

It seems feasible to solve the issue by averaging over many realisations of similar flow structures each averaged over short time (therefore noisy). This concept was already used by other authors (e.g., Duvall et al. 2006). As an example we note the possibility to measure the flows in many individual supergranules and to average over this sample in order to obtain the typical flow structure in an “average supergranule”, similar to Duvall & Birch (2010). The predicted error in the results after averaging scales as $\mathrm{1}\mathit{/}\sqrt{\mathrm{𝒩}}$, where is the size of the ensemble of independent flow realisations. This scaling law allows us to relax the noise constraint on the inversion and therefore regularise more strongly about the misfit and cross-talk terms. See Table 2 for required random error of inversions results.

For instance, let us assume that we average over 10⁴ independent inverted flows, each obtained by considering travel-time maps averaged over 6 h. Individual flow maps are noisy and therefore contain little useful information about underlying flows. This selection allows us to relax the constraint on the estimated error level by a factor of a 100 and therefore obtain a much better fit to the target function (see Fig. 27). With this set-up, the validation displayed in Fig. 15 shows almost perfect correspondence between expected v^tgt and inverted $\begin{array}{}{}_{\mathbf{˜}}\\ {v}\end{array}{}^{\mathrm{inv}}$ flows in the top 5.5 Mm of the convection zone. Inversion for v_z at 5.5 Mm with a FWHM of s_h = 15 Mm is an exception, because the selected target function cannot be matched by the averaging kernel. Different inversions for v_z at this depth are possible with greater averaging, e.g., with s_h = 25 Mm. It is still crucial to minimise cross-talk in the case of v_z inversions to retrieve the correct answer. Power spectra of inversion components at all discussed depths are displayed in Fig. 28.

7. Conclusions

We improved and validated an inversion algorithm based on a SOLA inversion approach. The formalism, algorithm and the code is universal and can, in principle, be used to invert for any quantity describing inhomogeneities in the solar plasma, provided that the corresponding sensitivity kernels are available. The code is also ready to be used for application to real measurements. It will become part of the helioseismic pipeline running at German Data Center for SDO. We plan to use it to analyse all available SOHO/MDI and SDO/HMI data in order to routinely provide a tomographic image of the structure of the solar upper convection zone.

The code is absolutely scalable allowing to include more independent measurements (and therefore more sensitivity kernels) in order to further refine the precision of the results. Thanks to the decoupling, the problem will still be solvable even using nowadays computers.

We improved the inversions by introducing additional terms that allow to control and minimise some sources of bias in the results. Most importantly, we minimise the cross-talk between individual flow components, which is crucial especially for v_z inversions. In principle, the formalism allows to minimise cross-talk between any selected quantities during the general inversion. The validation performed here silently assumes that the sensitivity kernels and noise covariance matrices are perfect. For the validation of the method and the code it is important if the sensitivity kernels and noise covariance matrices are solar-like. This allows us to study different sources of biases we may expect in the Sun. It is also very important that the sensitivity kernels used in the inversion contain all the details of the travel-time measurements, including the instrumental function affecting the solar oscillation power spectrum.

We found that by considering f to p₄ frequency-averaged modes and supergranular spatial scales it is possible to perform reliable and trustworthy flow inversions of the travel-time maps averaged over a few days in the top 3.5 Mm layer of the convection zone in the case of horizontal v_x and v_y components, and in the top 1 Mm in the case of vertical v_z component. Based on our experiment, we expect that using travel-time maps averaged over 4 days, it is possible to measure 3D velocities as weak as 10 m s^-1 at the surface and horizontal velocities having amplitude 20 m s^-1 in the top 5 Mm. We estimate, that using the travel-times averaged over 1 day, it should still be possible to measure all components of the supergranular flow at the surface and its horizontal components in the top 5 Mm. The cross-talk minimisation is crucial in order to measure the correct vertical velocity. Its presence in the results may explain the opposite sign of the vertical flow inversion discovered by Zhao et al. (2007).

By considering many flow realisations and statistical averaging we might go deeper in the convection zone to learn about the horizontal flow components and perform vertical flow inversions for the top few Mm depths.

Online material

Acknowledgments

This study was supported by the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007–2013)/ERC grant agreement #210949, “Seismic Imaging of the Solar Interior”, to PI L. Gizon (Milestone #5). Authors would like to thank J. Jackiewicz for providing us with the noise covariance matrices. The flow sensitivity kernels were computed using the code written by A. C. Birch deployed within the HELAS project at http://www.mps.mpg.de/projects/seismo/NA4. M.Š. acknowledges a partial support through the Grant Agency of Academy of Sciences of the Czech Republic under grant IAA30030808.

References

All Tables

All Figures