Recovery of subsurface profiles of supergranular flows via iterative inversion of synthetic travel times
^{1} Tata Institute of Fundamental Research, 400005 Mumbai, India
email: jishnubhattachrya@gmail.com
^{2} Max Planck Institute for Solar System Research, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
^{3} Institut für Astrophysik, GeorgAugustUniversität Göttingen, 37077 Göttingen, Germany
^{4} Center for Space Science, NYUAD Institute, New York University Abu Dhabi, 129188 Abu Dhabi, UAE
Received: 3 May 2017
Accepted: 28 July 2017
Aims. We develop a helioseismic inversion algorithm that can be used to recover subsurface vertical profiles of twodimensional supergranular flows from surface measurements of synthetic wave travel times.
Methods. We carried out seismic wavepropagation simulations with a twodimensional section of a flow profile that resembles an average supergranule and a starting model that only has flows at the surface. We assumed that the wave measurements are entirely without realization noise for the purpose of our test. We expanded the vertical profile of the supergranule stream function on a basis of Bsplines. We iteratively updated the Bspline coefficients of the supergranule model to reduce the traveltime differences observed between the two simulations. We performed the exercise for four different vertical profiles peaking at different depths below the solar surface.
Results. We are able to accurately recover depth profiles of four supergranule models at depths up to 8−10 Mm below the solar surface using f−p_{4} modes under the assumption that there is no realization noise. We are able to obtain the peak depth and the depth of the return flow for each model.
Conclusions. A basisresolved inversion performs significantly better than an inversion in which the flow field is inverted at each point in the radial grid. This is an encouraging result and might act as a guide in developing more realistic inversion strategies that can be applied to supergranular flows in the Sun.
Key words: Sun: helioseismology / waves / methods: numerical
© ESO, 2017
1. Introduction
Convective flows on the solar surface exhibit several length scales (Nordlund et al. 2009). Smallscale granules (~1 Mm) are well studied and characterized. The subsurface profile and physics behind relatively larger scale flows have eluded a convincing explanation. Supergranules, which are believed to be overturning flows spanning 35 Mm horizontally on average (Hathaway et al. 2000; Rieutord et al. 2008), remain one of the areas in which subsurface imaging has had limited success. Sketching a complete picture of a supergranule involves understanding the nature and magnitude of upflows near the cell center, downflows at the edges of the cell, horizontally diverging flows from the center of the cell toward the edge, and deeper return flows that might be present. There are several techniques used to study subsurface flows (see Gizon et al. 2010, and references therein), and from these we focus on timedistance seismology (Duvall et al. 1993). This approach allows us to relate traveltime shifts of seismic waves in the sun to surface and subsurface velocities, thereby setting up an inverse problem in which measurements of wave travel times can be used to estimate flow fields in the solar interior. Timedistance seismology has been widely used to recover flows in the Sun (Duvall & Gizon 2000; Zhao & Kosovichev 2003; Zhao 2004; Jackiewicz et al. 2008; Duvall & Hanasoge 2013; Švanda 2012), however when applied specifically to supergranules, the subsurface flow profiles have been hard to pin down.
Despite successful measurements of supergranular flows on the solar surface (Rieutord et al. 2008; Duvall & Birch 2010; Švanda et al. 2013), seismic studies have not been successful at consistently reproducing the subsurface profiles of supergranules. Duvall (1998) used correlations between inverted surface and deeper flows to determine the lower bound of the supergranule pattern where the convective cell overturns and obtained a depth of 8 Mm. Zhao & Kosovichev (2003) used timedistance seismology to invert MDI data and inferred that the depth of a supergranule was 15 Mm. Braun et al. (2004) applied phasesensitive holography to MDI data and concluded that detection of return flow below 10 Mm would prove to be a significant challenge because of contamination from neighboring supergranules. Woodard (2007) used Fourierspace correlations in the observed wave field and found power extending down to 6 Mm below the surface before noise took over; similar results were also obtained by Braun et al. (2007) via helioseismic holography and Jackiewicz et al. (2008) via timedistance seismology. Using measurements of horizontal flow divergences and vertical velocity obtained from the Solar Optical Telescope (SOT) on board the Hinode satellite, Rieutord et al. (2010) estimated the vertical scale height of supergranules to be 1 Mm. Duvall & Hanasoge (2013) used a raytheoretic forward modeling approach to fit centerannulus traveltime differences obtained with Gaussian models of vertical velocity to those obtained from a kinematic model of an “averaged supergranule” derived using Dopplergrams from the Helioseismic and Magnetic Imager (HMI: Schou et al. 2012) on board the Solar Dynamics Observatory (SDO) spacecraft. These authors estimated that the depth corresponding to peak vertical velocity to be 2.3 ± 0.9 Mm, where the value signified by ± represents the width of the model. Additional evidence for a shallow supergranule was obtained by Duvall et al. (2014).
The prevalence of disparate values and the dependence of results on the use of a specific technique calls for validation tests of seismic inversion algorithms. Dombroski et al. (2013) tested regularized leastsquare inversions via helioseismic holography measurements for a supergranulationlike flow, but found that the inferred vertical flow has significant errors throughout the computational domain. Švanda et al. (2011) used subtractive optimally localized averages (SOLA; Jackiewicz et al. 2008) and were able to recover threedimensional velocity fields from forwardmodeled, traveltime maps generated from simulations of solarlike convective flows. However, Švanda (2015) showed that reconstructed velocity fields do not produce wave travel times that match the observed travel times. DeGrave et al. (2014) tried validating timedistance SOLA inversions with realistic solar simulations and found that they could recover horizontal flows down to 5 Mm below the solar surface; but, unlike Švanda et al. (2011), DeGrave et al. were unable to infer vertical flows accurately. The authors attributed this to differences in measurement and analysis techniques. Hanasoge (2014) and Bhattacharya & Hanasoge (2016) tried a different appraoch, which used fullwaveform inversion (Tromp et al. 2010) to update iteratively a Cartesian twodimensional flow profile to minimize the traveltime misfit computed with respect to a model similar to the average supergranule from Duvall & Birch (2010). These authors found that seismic waves in their simulations were primarily sensitive to flow updates close to the solar surface and inversions focused on updating these layers at the expense of deeper layers. This negative result was surprising, especially for a noisefree inversion, in light of the previous studies by Švanda et al. (2011) and DeGrave et al. (2014), and urged a deeper probe into the reasons behind this mismatch.
In this work, we follow an approach similar to Bhattacharya & Hanasoge (2016), but pose the problem differently to avoid the interplay between the large number of parameters being determined and reduce the question to the fundamental one of seismic sensitivity to flows. We approach the question by considering the pedagogic exercise of inverting a twodimensional section of the averaged supergranule of Duvall & Hanasoge (2013), using seismic waves that are excited by sources located at specific spatial locations. While the setup is not directly comparable with solar observations, this serves as a computationally efficient starting point to validate fullwaveform inversion applied to the Sun. We project the supergranular flow model in a Bspline basis and solve an optimization problem to obtain the spline coefficients. We show that we are able to recover accurately the vertical profile of the averaged supergranule down to 8−10 Mm below the surface. This result might help to guide the construction of improved inversion strategies to study the subsurface profile of an averaged supergranule in the Sun.
2. Supergranule model
We considered kinematic models of temporally stationary supergranules in Cartesian coordinates. The entire analysis is twodimensional primarily for computational ease, although it provides us with an extra simplification that would be absent in three dimensions – that of a unidirectional stream function. We chose coordinates x = ^{(}x,z^{)}, where z is a vertical coordinate that increases in the direction opposite to gravity, and x denotes a horizontal direction with periodic boundary conditions. We set z = 0 at the solar surface, therefore negative values of z indicate depths below and positive values indicate heights above it.
A model of a supergranule can be described by its velocity field that is embedded in a steady solar background characterized by a onedimensional density profile , soundspeed profile , pressure , and acceleration due to gravity . The velocity profile of the supergranule is chosen to resemble a section through the averaged supergranule of Duvall & Hanasoge (2013); the differences arise because the computation is in Cartesian coordinates rather than cylindrical coordinates.
We enforced mass conservation (1)and derived the velocity field from a stream function as (2)Our model for the supergranule stream function is (3)where J_{1} is the Bessel function of order 1. The supergranule stream function is zero at the cell center, and peaks at a certain distance away from the center before falling to zero and reversing sign; the reversal in sign indicates a transition in the vertical velocity from upflows to downflows. The model is highly simplified compared to supergranules as observed on the solar surface because it ignores the impact of magnetic fields and other observed anomalous characteristics, such as a wavelike nature associated with supergranules (Gizon et al. 2003) and eastwest traveltime asymmetries (Langfellner et al. 2015). We fixed the horizontal length scales to R = 15 Mm and k = 2π/^{(}30 Mm^{)} and used several combinations of parameters to characterize the vertical profile. These parameters are listed in Table 1.
The flow field that we obtained from this stream function is We list the magnitude of the peak and surface velocities for these flow fields in Table 2. In subsequent analyses, we shall refer to this velocity field with the superscript “true”, i.e., as v^{true} and similarly for its components, to distinguish it from the flow velocity in the iteratively updated flow model. We apply superscripts “true” and “iter” to other parameters wherever necessary, indicating which model they correspond to.
Streamfunction parameters.
3. Inversion setup
Waves in the Sun are driven near the solar surface by turbulent convection associated with granules with most of the excitation taking place within 500 km of the photosphere (Stein & Nordlund 2001). Once generated, seismic waves propagate under the restoring forces applied by fluid pressure gradients and gravity. The wave displacement ξ(x,t) evolves according to the equation (6)where S represents sources that are exciting waves in the Sun. In our simulation, we chose eight sources located at different horizontal positions at a depth of 150 km below the surface. The sources represent “master pixels” (Tromp et al. 2010; Hanasoge et al. 2011), that is their locations are chosen so that the emanating waves sample the supergranule adequately. Each source fires independently of the others and produces waves that illuminate slightly different regions in the Sun. For both the true supergranule and the iterated supergranule, therefore, we computed eight different simulations in parallel, each of which ran for 4 h in solar time. The simulation box spans 800 Mm horizontally over 512 pixels, extends from 137 Mm below the surface to 1.18 Mm above it vertically, and is resolved with 300 pixels spaced uniformly in acoustic distance. We placed perfectly matched layers along the vertical boundaries to absorb waves effectively.
Fig. 1 Left panel: wave spectrum overlaid with the function used to filter waves corresponding to the radial order p_{2}. Middle and right panels: traveltime shifts between waves in the starting and true flow models for supergranule SG2, computed at all receivers for different radial orders. The horizontal location of the source is indicated by the dashed vertical line. The horizontal flow field of the supergranule is indicated by the colored patch in the background, where red indicates outflows away from the cell center and blue indicates inflows toward the cell center. 

Open with DEXTER 
Peak velocities, surface velocities, and depths for the four models considered.
We used the seismic wave propagation code SPARC (Hanasoge & Duvall 2007) to solve the wave equation in a convectively stabilized version of Model S (ChristensenDalsgaard et al. 1996). The model is onedimensional, satisfies hydrostatic balance, and is stabilized by patching an isothermal layer to Model S above 0.98 R_{⊙} (Hanasoge et al. 2006). The code SPARC computes seismic wave fields by solving Eq. (6) in the time domain using a fivestage, RungeKutta timestepping scheme with low dispersion and low dissipation (Hu et al. 1996). We computed the spatial derivatives with a sixthorder, compact finitedifference scheme (Lele 1992) in the vertical direction and Fourier decomposition in the horizontal direction.
We applied ridge filters to study the propagation of wave packets corresponding to individual radial orders. The ideal filtering technique has been a subject of some debate. Švanda (2013) found that inversions using a ridgefiltered approach produces results that are consistent with those using a phasespeed filtered approach, while DeGrave et al. (2014) found that inversions using ridgefiltered travel times do not compare favorably with phasespeedfiltered travel times; however we do not address this issue in the present work. Following Jackiewicz et al. (2008), we filtered the data by multiplying the wave spectrum by a function of the form F_{n}^{(}ν,k^{)}, where ν represents temporal frequency, k represents spatial frequency, and n represents the radial order. We constructed the filter function by separating modes corresponding to different radial orders with fourthorder polynomials of k. For each radial order, the spectral area enclosed between two such polynomials and is entirely included. We listed the polynomials used for each radial order in Table 3. The selected temporalfrequency band at each pixel in k is terminated with a quarter of a cosine function over two pixels on both the high and low edges to ensure a smooth fall off. Additionally low temporal frequency modes below 1.1 mHz are filtered out to remove contributions from weak gmodes that arise as an artifact of the convectively stabilized background.
We chose a group of pixels 200 km above the surface and marked them as receivers. We listed the horizontal locations of receivers for various radial orders in Table 3. The wide range of receiver locations combines both short and largedistance measurements, thereby utilizing waves that probe various depths beneath the solar surface. We recorded filtered waveforms with time for f, p_{1}, p_{2}, and p_{3} ridges at each receiver pixel for each simulation, for the p_{4} ridge, and for the case of SG4. An example of a spectrum, along with a filter to extract waves corresponding to the radial order p_{2} and measured travel times at receivers for the model SG2, is depicted in Fig. 1. We define the traveltime misfit (7)where refers to the waves emanating from the source s whose ridgefiltered travel time is measured at the rth receiver in presence of the true supergranule, is the travel time measured in presence of the iterated model for the same sourcereceiver locations and the same filter, and the sum extends over sourcereceiver pairs as well as different radial orders. We computed travel times following Gizon & Birch (2002) and we describe the technique in detail in Appendix A.
Nonlinear iterative timedistance inversions, as formulated in the context of helioseismology by Hanasoge (2014), revolve around reducing the misfit defined in Eq. (7) by sequentially improving a model of the supergranule. The scheme proceeds by relating the traveltime misfit to an update in the supergranule stream function through an integral relation as (8)where is the kernel whose value at any spatial point indicates sensitivity of wave travel times to local updates in the supergranule model. We used the adjoint source technique (Hanasoge et al. 2011) to compute the finitefrequency kernel . The steps involved in computing this kernel have been detailed in Hanasoge (2014) and Bhattacharya & Hanasoge (2016).
Sourcereceiver distances.
3.1. Basisresolved inversion
Previous attempts at inversions by Hanasoge (2014) and Bhattacharya & Hanasoge (2016) have focused on solving for the stream function at each spatial location. The number of parameters being determined was equal to the number of spatial grid points, i.e., we would have had 153 600 parameters for a 512 × 300 grid. We wondered whether the large size of the parameter space kept previous attempts from converging to the correct model. To answer this question, we posed the inverse problem differently and made the following assumptions:
 1.
The supergranule stream function is separable in x and z. This assumption is made keeping in mind that we are primarily interested in the depth of supergranules.
 2.
The value of the stream function at the surface and above is known. The justification is that the layers at and above the surface are directly observed and hence the flow velocities are measured (Gizon et al. 2000; Rieutord et al. 2008; Duvall & Birch 2010).
We reduced the parameter space further by expanding the vertical profile of ψ^{true}(x) in a basis of Bsplines and inverting for the coefficients close to the surface. We expanded the vertical profile in a basis of quadratic Bsplines with a given set of knots {t} as (10)where β_{i} represent the Bspline coefficients and k = 2 indicates quadratic splines. The Bsplines are ordered such that the index i = 0 corresponds to the Bspline function that peaks the deepest, while the index i = N−1 corresponds to the index that peaks close to the upper boundary of our computational domain. The approximate equality is to be understood as the best fit in a leastsquares sense, since we chose a set of smoothing splines instead of interpolating splines. We describe the spline expansion in detail in Appendix B. Such an approach is only applicable to an ensemble averaged model of a supergranule, where the flow profile is expected to be smooth and representable with relatively few splines.
We split the coefficients into two groups: those above the surface and those below the surface. Assuming that the coefficient with index m corresponds to the Bspline function that peaks at the solar surface, we rewrote Eq. (10) as (11)We chose the surface profile and used Eq. (9) to obtain the starting model of our supergranule, i.e., (12)The starting flow profile in our inversion is the same as the true flow above the surface and falls to zero continuously just below the surface. We can represent the vertical profile of the starting supergranule stream function on a basis of splines by setting the coefficients of Bsplines below the surface to zero, i.e., We iteratively updated this model. Therefore, at the first step of the inversion we set ψ^{iter} = ψ^{start}, or equivalently . This choice is different from that made by Hanasoge (2014) and Bhattacharya & Hanasoge (2016), where the starting model had no flows. We used the same set of knots { t } to represent the inverted model as those that we had used to expand g^{true}(z) in Eq. (10). The inversion is carried out to obtain the spline coefficients that lie below the solar surface.
We substituted the spline expansion of the iterated stream function in Eq. (8) to obtain kernels in spline space as (15)The discrete kernels K_{i} indicate the sensitivity of travel times to individual Bspline coefficients. We used the BroydenFletcherGoldfarbShanno algorithm (BFGS, Nocedal & Wright 2006) to iteratively update our model of the supergranule flow profile and reduce the traveltime misfit.
Fig. 2 True and inverted flow velocity for SG2 from Table 1 with data and model misfits. The top left panel indicates true v_{x}, the top center indicates inverted v_{x}, and the top right indicates misfit in v_{x} as a function of depth (Eq. (18)). The middle left panel indicates true v_{z}, middle center indicates inverted v_{z}, middle right indicates misfit in v_{z} as a function of depth (Eq. (19)). The bottom left panel indicates data misfit from Eq. (7); the bottom center indicates model misfit from Eq. (20) for the stream function ψ and the two components of velocity. 

Open with DEXTER 
Fig. 3 True and inverted profiles for the four supergranule Gaussian profiles from Table 1. Each row corresponds to one model, where the left panel represents the Bspline coefficients, the middle panel shows the vertical profile of the stream function, and the rightmost panel depicts the vertical profile of the horizontal component of the flow velocity. The profile of vertical flow is not plotted, but is similar to the stream function. In the leftmost panel for each row, bars represent Bspline coefficients for the true stream function; white squares represent spline coefficients above the surface, which are clamped to the value in the true model; and black circles represent coefficients for the inverted solution. In the middle and right panels of each row, gray solid lines represents vertical profiles of the true models and black circles denote the profile for the inversion result. 

Open with DEXTER 
Fig. 4 Left: peak depths of the inverted stream function vs. those of the true model for each of the different supergranules from Table 1. The dotted line indicates the ideal recovered peak depth, i.e., that of the true model. Right: reversal depths of v_{x} in the inverted models against those in the true models. 

Open with DEXTER 
3.2. Regularization
The simulated spectrum (left panel in Fig. 1) features discernible modal ridges from f to p_{6}, out of which we used ridges up to p_{4} for our study. Restricting ourselves to a small set of modes imposes a limit on the depth up to which we can infer flows, beyond this waves have limited sensitivity to flows. We therefore computed the knots required for the basis expansion in Eq. (10) and impose a lower cutoff. The depth of the cutoff is governed by the modes used in the inversion for each model, but this depth is chosen to be deep enough to ensure that the entire flow profile is contained within the spatial range.
We can estimate the depths that seismic waves probe by computing the asymptotic lower turning point (Giles 2000). The turning point for a wave that has the maximum power in the p_{3} ridge, corresponding to k_{x}R_{⊙} = 533 and temporal frequency ν = 4.5 mHz in the simulation, is 9 Mm. We therefore expected inferences using radial orders f−p_{3} to be accurate down to this depth. This also indicates that waves from the p_{4} and higher radial orders might be necessary to infer flows deeper down. The value of the lower cutoff and number of spline coefficients used for each model are listed in Table B.1. We ensured that the iterated flow model falls smoothly to zero at the lower cutoff by multiplying the ith Bspline coefficient by a factor of 1 / (1 + exp(−(i−i_{2}) / 0.2)), where i_{2} is the index of the coefficient that lies 2 Mm above the lower cutoff. For models SG1, SG2, and SG3 we obtained i_{2} = 1, indicating that the two deepest coefficients (i = 0and i = 1) are suppressed by factors of 150 and 2, respectively, whereas for SG4 we obtain i_{2} = 0, indicating that the value of the deepest coefficient (i = 0) is reduced by a factor of 2.
Previous analysis by Bhattacharya & Hanasoge (2016) had used spatial smoothing to reduce high spatialfrequency variation in the numerically computed sensitivity kernel. This is not strictly necessary in our approach and we found minor differences by including smoothing.
4. Results and discussion
The inversion in our analysis is a series of forward simulations, followed by optimization in the parameter space of the flow. Each stage in the iterative optimization proceeds by reducing the traveltime misfit in Eq. (7). We plot the travel time for different radial orders for the model SG2 in Fig. 1. We compare the traveltime shift with the horizontal flow that is indicated by the colored patch. The traveltime shift is a measure of how much the waveform in the starting model is delayed with respect to that in the true model, negative values indicating that the wavepacket in the starting model arrives earlier at a receiver in relation to the true model.
Each model is updated by iteratively reducing the traveltime misfit using Eq. (8). In solar traveltime measurements, error bars arising from realization noise provide a natural stopping point for iterations; in the absence of noise we iterate until the relative change in travelmisfit falls below 0.1%. We quantify the efficacy of the inversion by defining model misfits for the stream function and the components of the flow. The flow velocities are related to derivatives of stream function through Eq. (2), so the misfit in components of flow, when computed at each depth, differs from that for the stream function. The misfit κ for each parameter is defined as the normalized square of the difference between true and iterated models evaluated as a function of depth by averaging over the horizontal direction x, i.e., for the stream function ψ we obtain (16)where the horizontal integral in the denominator is evaluated at the depth at which the true model reaches its peak. This ensures that the maximum misfit for the starting model is normalized to one. The notation here indicates that the misfit κ is computed for the parameter in square brackets and is evaluated as a function of vertical layer z. Using the separability condition in Eq. (9), we can express κ[ψ] in terms of the vertical profile as (17)We define analogous misfit functions for the two components of flow velocity, where the expressions are Alongside studying misfit as a function of depth, we consider the normalized L_{2} norm of differences between true and iterated models integrated over the entire space, defined as (20)to gain insight into the degree of improvement to the model after each iteration. We compare the inverted flow velocity field with the profile of the model SG2 in Fig. 2 and analyze the model misfit. We find that the inversion result matches the true model reasonably well, where the stream function misfit κ_{L2}[ψ] is around 0.01% after the final iteration.
We plot the inverted vertical profiles of all the models in Fig. 3. One question we ask in this paper is whether seismic waves can estimate the depth of supergranules. We plot the expected and inferred peak depths for each model in Fig. 4. We find that we recover the peak depths accurately for all the models. It is encouraging that we are able to extract Gaussian profiles up to a depth of 8−10 Mm, that is beyond the 6 Mm limit found by previous seismic inferences in the presence of realization noise (Braun et al. 2007; Woodard 2007). We were however unable to retrieve the profile for a deeper model with a peak depth of 14 Mm with f−p_{4} modes. It might be interesting to see if the introduction of higher p modes makes a difference.
Aside from peak depth, another important parameter that is used to estimate the depth of supergranules is the layer at which the horizontal flow reverses direction. All the models that we study have a reversal in v_{x}, in keeping with the assumption of a steady convective cell. We compare the true and inferred reversal depths in Fig. 4. We find that the inversion reproduces comparable values, although the exact profile of the return flow is not accurately captured. Whether the reversal actually takes place is subject to debate, since it is not detected in the study by Woodard (2007) and is suggested to be spurious by Švanda (2013), DeGrave et al. (2014). Our results seem to indicate that the depth would be captured correctly for shallow models if the flow reverses provided that the inference is not limited by noise.
In this work, we used iterative forward modeling to obtain the bestfit flow model given traveltime measurements at the solar surface. This approach is inherently nonlinear, as it requires recomputing the sensitivity kernel after each iteration. From Fig. 2, we see that a substantial drop in model misfit takes place in the first few iterations. This makes it interesting to compare this approach with a linear inversion; given the small parameter space, it might be possible to pinpoint the differences in the inverted flow arising from the two approaches. We also have not considered any noise associated with traveltime measurements, so this analysis cannot be directly applied to solar measurements. It would be interesting to study the extent to which the inferences are affected in presence of noise and develop a slightly modified technique that accounts for realistic nose covariance matrices (Gizon & Birch 2004; Švanda et al. 2011). It was also pointed out by DeGrave et al. (2014) that validation tests with frozen nonmagnetic flow fields might be too idealized a scenario compared to studying a Doppler time series obtained from the Sun. Therefore an extension of this work to flows present in realistic solar simulations might be in order.
Acknowledgments
S.M.H. acknowledges support from Ramanujan fellowship SB/S2/RJN73/2013, the MaxPlanck partner group program, and thanks the Center for Space Science, New York University at Abu Dhabi. J.B. acknowledges the financial support provided by the Department of Atomic Energy, India.
References
 Bhattacharya, J., & Hanasoge, S. M. 2016, ApJ, 826, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Braun, D. C., Birch, A. C., & Lindsey, C. 2004, in SOHO 14 Helio and Asteroseismology: Towards a Golden Future, ed. D. Danesy, ESA Spec. Publ., 559, 337 [Google Scholar]
 Braun, D. C., Birch, A. C., Benson, D., Stein, R. F., & Nordlund, Å. 2007, ApJ, 669, 1395 [NASA ADS] [CrossRef] [Google Scholar]
 ChristensenDalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 DeGrave, K., Jackiewicz, J., & Rempel, M. 2014, ApJ, 788, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Dierckx, P. 1993, Curve and Surface Fitting with Splines (New York, NY, USA: Oxford University Press, Inc.) [Google Scholar]
 Dombroski, D. E., Birch, A. C., Braun, D. C., & Hanasoge, S. M. 2013, Sol. Phys., 282, 361 [NASA ADS] [CrossRef] [Google Scholar]
 Duvall, Jr., T. L. 1998, in Structure and Dynamics of the Interior of the Sun and Sunlike Stars, ed. S. Korzennik, ESA Spec. Publ., 418, 581 [Google Scholar]
 Duvall, Jr., T. L., & Birch, A. C. 2010, ApJ, 725, L47 [NASA ADS] [CrossRef] [Google Scholar]
 Duvall, Jr., T. L., & Gizon, L. 2000, Sol. Phys., 192, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Duvall, Jr., T. L., & Hanasoge, S. 2013, Sol. Phys., 287, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Duvall, Jr., T. L., Jefferies, S., Harvey, J., & Pomerantz, M. 1993, Nature, 362, 430 [NASA ADS] [CrossRef] [Google Scholar]
 Duvall, Jr., T. L., Hanasoge, S. M., & Chakraborty, S. 2014, Sol. Phys., 289, 3421 [NASA ADS] [CrossRef] [Google Scholar]
 Giles, P. M. 2000, Ph.D. Thesis, Stanford University [Google Scholar]
 Gizon, L., & Birch, A. C. 2002, ApJ, 571, 966 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., & Birch, A. C. 2004, ApJ, 614, 472 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., Duvall, Jr., T. L., & Larsen, R. M. 2000, JA&A, 21, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Gizon, L., Duvall, T. L., & Schou, J. 2003, Nature, 421, 43 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Gizon, L., Birch, A. C., & Spruit, H. C. 2010, ARA&A, 48, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M. 2014, ApJ, 797, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M., & Duvall, Jr., T. L. 2007, Astron. Nachr., 328, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M., Larsen, R. M., Duvall, Jr., T. L., et al. 2006, ApJ, 648, 1268 [NASA ADS] [CrossRef] [Google Scholar]
 Hanasoge, S. M., Birch, A., Gizon, L., & Tromp, J. 2011, ApJ, 738, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Hathaway, D. H., Beck, J. G., Bogart, R. S., et al. 2000, Sol. Phys., 193, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, F. Q., Hussaini, M. Y., & Manthey, J. L. 1996, J. Comput. Phys., 124, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Jackiewicz, J., Gizon, L., Birch, A. C., & Duvall, Jr., T. L. 2007, ApJ, 671, 1051 [NASA ADS] [CrossRef] [Google Scholar]
 Jackiewicz, J., Gizon, L., & Birch, A. C. 2008, Sol. Phys., 251, 381 [NASA ADS] [CrossRef] [Google Scholar]
 Langfellner, J., Gizon, L., & Birch, A. C. 2015, A&A, 579, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lele, S. K. 1992, J. Comput. Phys., 103, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Nocedal, J., & Wright, S. J. 2006, Numerical Optimization, 2nd edn. (New York: Springer) [Google Scholar]
 Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Liv. Rev. Sol. Phys., 6, 2 [Google Scholar]
 Rieutord, M., Meunier, N., Roudier, T., et al. 2008, A&A, 479, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M., Roudier, T., Rincon, F., et al. 2010, A&A, 512, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Sol. Phys., 275, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Stein, R. F., & Nordlund, Å. 2001, ApJ, 546, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Švanda, M. 2012, ApJ, 759, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Švanda, M. 2013, ApJ, 775, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Švanda, M. 2015, A&A, 575, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Švanda, M., Gizon, L., Hanasoge, S. M., & Ustyugov, S. D. 2011, A&A, 530, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Švanda, M., Roudier, T., Rieutord, M., Burston, R., & Gizon, L. 2013, ApJ, 771, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Tromp, J., Luo, Y., Hanasoge, S., & Peter, D. 2010, Geophys. J. Int., 183, 791 [NASA ADS] [CrossRef] [Google Scholar]
 Woodard, M. F. 2007, ApJ, 668, 1189 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, J. 2004, Ph.D. Thesis, Stanford University [Google Scholar]
 Zhao, J., & Kosovichev, A. G. 2003, in GONG+ 2002. Local and Global Helioseismology: the Present and Future, ed. H. SawayaLacoste, ESA Spec. Publ., 517, 417 [Google Scholar]
Appendix A: Traveltime measurements
We measured traveltime shifts by minimizing the squared difference between wave displacements in the true and starting supergranule simulations. Given wave displacements ξ^{true} and ξ^{iter} recorded at a receiver at x_{r} and filtered to obtain the wave packet corresponding to a specific radial order, we defined a misfit (A.1)where is a window function that encloses the wave packets and isolates them from artifacts that might arise from spatiotemporal periodicity assumed in the simulation. Specifically, we chose the functional form of to be a box function, that is one over the range of the wave packet and zero outside. We plot one example of a measured wave packet and window function in Fig. A.1 (left panel). The traveltime shift is defined as the value of τ that minimizes (Gizon & Birch 2002, 2004). For bandlimited waveforms sampled beyond twice their Nyquist frequency, this can be computed by equating the temporal derivative of to zero and solving for τ. Expanding η in terms of the timeshift τ, suppressing the explicit dependence on the coordinates x_{r} and t, representing the order of derivative using superscripts in parentheses and referring to as δξ, we obtained \newpage
Retaining terms until quadratic order and solving leads to a traveltime shift given by (A.2)Equation (A.2) is in the form of traveltime shift defined by Gizon & Birch (2002), where the shift δτ is linear in the displacement difference δξ. The linear dependence of traveltime shifts on δξ is important for consistent computation of sensitivity kernels in the first Born approximation. At high flow velocities, however, it is possible that this linear relationship fails to remain a good approximation (Jackiewicz et al. 2007; DeGrave et al. 2014). For seismic waves with frequency ω, the deviation from linearity is noticeable if the measured traveltime shift δτ satisfies ωδτ ≪ 1. We carried out validation tests for by artificially timeshifting a simulated wave field by a typical value measured at surface for each supergranule model, and trying to recover the shift from the series expansion of truncated at various orders in τ. We plot the result of one such validation test in Fig. A.1 (right panel). We found that travel times computed using the linear approximation for the various supergranule models are within 5% of the expected value; the accuracy of the estimate improves with an increase in degree of the truncation. This is the error in traveltime measurement at the first iteration in our inversion. Subsequent iterations improve the flow model and lead to a significant decrease in δτ, consequently the error in measuring traveltime shifts is also reduced. This ensures that the inverted flow is not affected by the error in measured travel time for the starting model. It might be necessary to take this error into consideration for a linear inversion.
Fig. A.1 Left: one example of an fmode wave packet measured at a receiver located at x_{r} = 30 Mm for the supergranule SG1. The extent of the window function is denoted by vertical dashed lines. Right: relative error in estimated traveltime shifts using wave packets corresponding to radial orders f and p_{1} (indicated by different colors) as a function of different degree of truncation of Eq. (A.2) (indicated by different symbols), for different supergranule models. Truncating Eq. (A.2) to quadratic order results in a linear relation between traveltime shift δτ and wave displacement ξ^{iter} (plotted with diamonds). 

Open with DEXTER 
Appendix B: Spline expansion
We expanded the supergranule model from Eq. (3) in a basis of Bsplines following Eq. (10) to obtain a set of coefficients that we determined. The first step to this expansion was to obtain a set of knots that govern the Bsplines. The choice of knots might be important for the inversion to succeed, however further study needs to be carried out to establish this. We chose knots by applying the Dierckx algorithm (Dierckx 1993) to the true supergranule, but alternate choices of knots derived from an independent, vertically stratified parameter such as sound speed may also be used. We computed the knots and coefficients using the “scipy.interpolate” package of the python programming language, which internally calls the Fortran library FITPACK. The module computes the spline fit by evaluating an optimal set of knots and coefficients, ensuring that the squared L_{2} norm of the difference between the data being fit and the spline approximant falls below a specified smoothing factor. Reducing the value of this factor improves the fit; this is achieved by updating the set of knots followed by reevaluating the expansion coefficients. We carried out our inversion in a space spanned by the spline coefficients, therefore the specific choice of a approximant is a
trade off between the quality of fit and the number of parameters used to obtain the fit. This is why we chose smoothing splines over interpolating ones, since the latter involves a similar fit that is evaluated without any smoothing and produces a set of knots similar in size to the number of grid points, while the former can be tuned to significantly trim down the size of this set. We selected the smoothing factor through experimentation to obtain an accurate representation of the stream function in the Bspline basis while restraining the number of coefficients to around 10. We listed the smoothing parameters and knots for each of the supegranule models in Table B.1. We decided upon quadratic splines as a compromise between a cubicspline fit that is more oscillatory and a linearspline fit that is not as smooth and results in a larger parameter space.
We plot the functional form of the Bspline functions for model SG2 in Fig. B.1 (left panel). We plot the functional form of the supergranule from Eq. (10) and the spline approximation to it in Fig. B.1 (middle panel), and we plot the error in the approximation in the right panel. We found that smoothing splines are a reasonably good representation of this supergranule model beneath the solar surface.
Fig. B.1 Left: Bspline functions used in the expansion of the supergranule model SG2. The specific choice of knots might be important for the inversion to converge to the correct model. The density of knots with depth is reflective of the degree of stratification. The Bsplines are not normalized by the corresponding coefficients. Middle: supergranule model from Eq. (3) (gray solid line) and the smoothing spline approximation to this model (black circles). Right: error in the spline approximation with depth. 

Open with DEXTER 
Bspline parameters used in the inversion.
All Tables
All Figures
Fig. 1 Left panel: wave spectrum overlaid with the function used to filter waves corresponding to the radial order p_{2}. Middle and right panels: traveltime shifts between waves in the starting and true flow models for supergranule SG2, computed at all receivers for different radial orders. The horizontal location of the source is indicated by the dashed vertical line. The horizontal flow field of the supergranule is indicated by the colored patch in the background, where red indicates outflows away from the cell center and blue indicates inflows toward the cell center. 

Open with DEXTER  
In the text 
Fig. 2 True and inverted flow velocity for SG2 from Table 1 with data and model misfits. The top left panel indicates true v_{x}, the top center indicates inverted v_{x}, and the top right indicates misfit in v_{x} as a function of depth (Eq. (18)). The middle left panel indicates true v_{z}, middle center indicates inverted v_{z}, middle right indicates misfit in v_{z} as a function of depth (Eq. (19)). The bottom left panel indicates data misfit from Eq. (7); the bottom center indicates model misfit from Eq. (20) for the stream function ψ and the two components of velocity. 

Open with DEXTER  
In the text 
Fig. 3 True and inverted profiles for the four supergranule Gaussian profiles from Table 1. Each row corresponds to one model, where the left panel represents the Bspline coefficients, the middle panel shows the vertical profile of the stream function, and the rightmost panel depicts the vertical profile of the horizontal component of the flow velocity. The profile of vertical flow is not plotted, but is similar to the stream function. In the leftmost panel for each row, bars represent Bspline coefficients for the true stream function; white squares represent spline coefficients above the surface, which are clamped to the value in the true model; and black circles represent coefficients for the inverted solution. In the middle and right panels of each row, gray solid lines represents vertical profiles of the true models and black circles denote the profile for the inversion result. 

Open with DEXTER  
In the text 
Fig. 4 Left: peak depths of the inverted stream function vs. those of the true model for each of the different supergranules from Table 1. The dotted line indicates the ideal recovered peak depth, i.e., that of the true model. Right: reversal depths of v_{x} in the inverted models against those in the true models. 

Open with DEXTER  
In the text 
Fig. A.1 Left: one example of an fmode wave packet measured at a receiver located at x_{r} = 30 Mm for the supergranule SG1. The extent of the window function is denoted by vertical dashed lines. Right: relative error in estimated traveltime shifts using wave packets corresponding to radial orders f and p_{1} (indicated by different colors) as a function of different degree of truncation of Eq. (A.2) (indicated by different symbols), for different supergranule models. Truncating Eq. (A.2) to quadratic order results in a linear relation between traveltime shift δτ and wave displacement ξ^{iter} (plotted with diamonds). 

Open with DEXTER  
In the text 
Fig. B.1 Left: Bspline functions used in the expansion of the supergranule model SG2. The specific choice of knots might be important for the inversion to converge to the correct model. The density of knots with depth is reflective of the degree of stratification. The Bsplines are not normalized by the corresponding coefficients. Middle: supergranule model from Eq. (3) (gray solid line) and the smoothing spline approximation to this model (black circles). Right: error in the spline approximation with depth. 

Open with DEXTER  
In the text 