Testing Helioseismic-Holography Inversions for Supergranular Flows Using Synthetic Data

Supergranulation is one of the most visible length scales of solar convection and has been studied extensively by local helioseismology. We use synthetic data computed with the Seismic Propagation through Active Regions and Convection (SPARC) code to test regularized-least squares (RLS) inversions of helioseismic holography measurements for a supergranulation-like flow. The code simulates the acoustic wavefield by solving the linearized three-dimensional Euler equations in Cartesian geometry. We model a single supergranulation cell with a simple, axisymmetric, mass-conserving flow. The use of simulated data provides an opportunity for direct evaluation of the accuracy of measurement and inversion techniques. The RLS technique applied to helioseismic-holography measurements is generally successful in reproducing the structure of the horizontal flow field of the model supergranule cell. The errors are significant in horizontal-flow inversions near the top and bottom of the computational domain as well as in vertical-flow inversions throughout the domain. We show that the errors in the vertical velocity are due largely to cross talk from the horizontal velocity.


Introduction
Wave travel times are influenced by the properties of the medium in which they propagate; helioseismic holography Braun, 1990, 2000) can be used to measure travel-time shifts, and is therefore useful for probing the interior structure of the Sun Gizon, Birch, and Spruit, 2010).
Supergranulation is a visible pattern on the solar surface often interpreted as a scale of solar convection. Despite extensive research by local helioseismology, there remains substantial uncertainty regarding the characteristic scales and fundamental processes responsible for the creation and sustainment of supergranulation patterns (e.g. Gizon and Birch, 2005;. Supergranulation cells were first detected by Hart (1954;1956); observations and analysis of data over a period of five decades have revealed average cell sizes in the range 15 -30 Mm (e.g. Hart, 1956;Leighton, Noyes, and Simon, 1962;Duvall, 1980;Hathaway, 1992;Hagenaar, Schrijver, and Title, 1997;De Rosa and Toomre, 2004;Hirzberger et al., 2008), with discrepancies in estimated scale due largely to differences in analysis techniques. Similarly, direct inferences of the velocity field have produced a range of horizontal-flow estimates (e.g. Hart, 1954;Simon and Leighton, 1964;Hathaway et al., 2002). , motivated by observed scale-dependency in surface kinetic-energy distributions, performed a spectral analysis of velocities estimated using granule tracking and Doppler measurements. They measured a velocity of 300 m s −1 at a scale of 36 Mm. Estimating vertical velocity has proven more difficult due to the noise contribution and low signal strength, particularly near the boundaries of supergranulation cells . Results from SOHO/MDI data (Hathaway et al., 2002) and Hinode/SOT data  produced estimates of 30 m s −1 . Duvall and Birch (2010) estimated cell-center vertical flow to be about 10 m s −1 using direct Doppler measurements averaged over 1100 supergranules.
Several theories have been advanced to explain the development of supergranulation flow patterns at the observed preferential scale, although confirmation or rebuttal of these theories is largely lacking. An early theory proposed by Simon and Leighton (1964) attributed supergranular scales to convective instability from the recombination of ionized helium. More recently, Rieutord et al. (2000Rieutord et al. ( , 2001 proposed a model whereby exploding granules trigger large-scale instability of the granular flow, leading to supergranulation. Ploner, Solanki, and Gadun (2000) used a 1D model to similarly show that the interaction and merging of individual granular plumes may determine supergranular scales. Rast (2003) demonstrated that large spatial and long temporal supergranular scales arise in a simplified advective model that simulates the interaction of many small-scale and short-lived granular downflow plumes. Crouch, Charbonneau, and Thibault (2007) were able to simulate correlations between cell size and magnetic activity using a model based on a random-walk approximation to the dispersal and interaction of small-scale magnetic elements at the solar surface. Qualitatively, they found that this process produced supergranule-like spatial patterns.
Debate continues regarding the structure of supergranule flows below the solar surface (e.g. Duvall, 1998;Zhao and Kosovichev, 2003;Woodard, 2007;Sekii et al., 2007), largely due to the difficulties associated with estimating relevant quantities using helioseismology (e.g. Gizon and Birch, 2005;. These challenges are compounded by the sound speed and density stratification, which effectively decrease resolution and signal-to-noise ratio with depth, making the detection of a return flow elusive. Kosovichev and Duvall (1997) first used the time-distance helioseismology technique (Duvall et al., 1993) to make measurements of supergranular flow patterns from Doppler measurements. Zhao and Kosovichev (2003) also used time-distance helioseismology to infer supergranule flow fields; they found that horizontal flows could be derived reliably within a few megameters of the surface, but that vertical flows remained uncertain. They determined supergranulation to be a relatively deep phenomenon with a full pattern depth of up to 15 Mm and a return flow measured below a depth of 5 -6 Mm. Braun, Birch, and Lindsey (2004), using helioseismic holography, attributed the measured return flow to signal leakage, and concluded that supergranulation is a relatively shallow phenomenon. Woodard (2007) used a forward model to perform subsurface inversions of Doppler data, however detection was limited to 4 -5 Mm below the photosphere. Jackiewicz, Gizon, and Birch (2008) used a novel 2+1D optimally localized averaging (OLA) method to study supergranulation in the top few Mm. Braun et al. (2007) and Zhao et al. (2007) used synthetic data from numerical simulations of wave propagation through supergranule-like flows to test the performance of helioseismic techniques. The former compared model travel times with those computed from synthetic observations using helioseismic holography. The latter computed inversions using time-distance helioseismology and compared the results to the modeled flows. Both studies documented limitations in the ability of the techniques to detect the full extent of the flow fields. Svanda et al. (2011) tested the OLA inversion technique of Jackiewicz et al. (2012) on synthetic travel times and showed improved ability to infer vertical flow fields.
In this study, we simulate wave propagation through a simple kinematic model of a supergranule flow pattern and make helioseismic holography measurements from observations of the resulting velocity field. The aim is to test the quality of regularized least squares (RLS: e.g. Zhao, Kosovichev, and Duvall, 2001;Kosovichev, 1996, in the context of local helioseismology) inversions and gain insight into the limitations of this commonly used technique in local helioseismology by comparing the calculations with a known simulated flow field. We find that the RLS technique applied to helioseismic holography measurements is able to infer some general features of the horizontal field, but fails to infer the vertical field throughout the domain. Errors in vertical-flow inversions have previously been attributed to cross-talk effects (e.g. Zhao et al., 2007;Jackiewicz et al., 2012); herein we show the individual contributions from the divergent flow that comprise these effects.

Numerical Algorithm
We used the SPARC code (Hanasoge et al., 2008 to solve the three-dimensional linearized Euler equations in Cartesian geometry. The code computes derivatives in the vertical direction using sixth-order compact finite differences (Lele, 1992), while derivatives in the horizontal directions are computed spectrally with periodic boundaries. An absorbing sponge is used to damp wave reflection at the top and bottom boundaries .
Wave excitation is applied 200 km below the photosphere via a source term in the momentum equation (Hanasoge, Duvall, and Couvidat, 2007). The forcing is statistically uniform in the horizontal directions and Gaussian in the vertical direction with FWHM 200 km. In the spectral domain, the forcing distribution is localized to within the wavenumber range k = 0 − 2 rad Mm −1 and has a frequency maximum at 3 mHz.
Two different 24-hour (solar time) simulation cases were performed. The first case used the supergranule model with velocity vector v = v xx + v yŷ + v zẑ , where x and y define the horizontal coordinate system and z defines vertical distance from the photosphere. The second simulation case used no background flows, which allowed for application of the "noise subtraction" technique in the analysis (Werne, Birch, and Julien, 2004).

Supergranulation Model
The static background model is a convectively stabilized (CSM A: Schunker et al., 2011) variant of Model S (Christensen-Dalsgaard et al., 1996, which specifies the time-invariant components of the density [ρ], pressure [p], first adiabatic index [Γ 1 ], gravity [g], and sound speed [c s ]. The computational box is a 300 3 grid, spanning nonuniformly from z bot = −25 Mm to z top = 2.5 Mm in the vertical direction. The grid spans 100 Mm uniformly in the x-and y-directions. We model a single supergranule with a simple axisymmetric mass-conserving flow ( Figure 1) prescribed by the curl of a potential function where, for the sake of simplicity, we assume that the potential is separable as with functions h(z) and F (r) defining the vertical and radial dependance of the potential. The variables φ, r, and z denote the angular, radial, and vertical cylindrical coordinates, respectively. The formulation of Equation 1 is chosen so that mass conservation ∇ · ρv = 0 is automatically enforced. We chose the radial distribution F (r) to be of the form where J 1 is the first order Bessel function of the first kind, k = 2π/30 rad Mm −1 , and R = 15 Mm. This formulation forces decaying reversal of flow direction with increasing radial distance from the axial center. This functional form is similar to the flow associated with the "average supergranule" as seen by, e.g. Duvall and Hanasoge (2012). From Equation (1), the radial and vertical mass-flux distributions, respectively, are We chose to model the v r component of flow as a sum of outflow [v out ] and inflow [v in ] Gaussian distributions, such that where α 1 = 250 m s −1 , α 2 = 6.27 m s −1 , z 1 = 0.2 Mm, z 2 = −15 Mm, and D 1 = D 2 = 5 Mm. Using mass conservation, the coefficient α 2 is calculated from the known density distribution ρ and v out with chosen coefficient α 1 . The potential function A φ can then be calculated directly from In obtaining Equation (7), we set the integration constant equal to zero since we want the vertical velocity v z to be zero at the lower boundary (see Equation 8 below). With A φ fully prescribed, the vertical flow [v z ] is calculated from Equation (5): Note that due to the large density gradient with solar depth, outflow velocities near the photosphere are much larger than return flow velocities at greater depth in order to conserve mass flux ρv ( Figure 1).

Travel-time Measurements
Starting from the vertical velocity taken at a height of 200 km, we used helioseismic holography (Lindsey and Braun, 2000) to measure center-quadrant local-control correlations. These correlations are analagous to center-quadrant time-distance correlations that have traditionally been used to measure timedistance travel times (e.g. Gizon and Birch, 2005;Gizon, Birch, and Spruit, 2010). We used the quadrant geometry and ridge filters described by Braun and Birch (2008). In addition to the ridge filters, we also applied band-pass frequency filters to isolate particular ranges in frequency. The filters had central frequencies of 2.75, 3.00, 3.25, 3.50, 3.75, 4.00, 4.25, and 4.50 mHz and all had widths of 0.25 mHz. We use δτ x (r) (δτ y (r)) to denote x (y) direction center-quadrant travel-time shifts and δτ oi = δτ out (r) − δτ in (r) to denote out minus in centerquadrant travel-time shifts, where δτ out (r) (δτ in (r)) represent travel-time shifts of the observed outgoing (incoming) waves. Figure 2 shows the power spectrum of the vertical velocity at a height of 200 km above the photosphere in the simulation without any imposed flows. The resonant frequencies of the simulation are close to those from model S. As a result, we apply the ridge filters and holography Green's functions of Braun and Birch (2008) without modification. The travel-time measurements shown in Figures 3 and 4 were produced by subtracting the travel-time measurements from the simulation without flows (i.e., "noise subtraction"; Werne, Birch, and Julien, 2004). Figure 3 shows center-quadrant travel-time differences [δτ x ] for a range of radial orders and frequencies for the case where the synthetic data is produced from the supergranulation model. The flow produces travel-time shifts that generally decrease in amplitude with increasing radial order and increase in amplitude with increasing frequency (i.e., the travel-time shifts decrease with increasing horizontal phase speed). The spatial pattern is much the same in all cases and is what would be expected given the flow geometry in Figure 1. Figure 4 shows the δτ oi travel-time differences resulting from the supergranulation model. The general trend of the travel-time shifts with radial order and frequency band is consistent with that in Figure 3, although the δτ oi measurements feature a larger signal-to-noise ratio.

Flow Inversions
The inverse problem is to estimate the flow field defined by the supergranule model given the travel-time measurements in Section 3. Here we show some example inversions in which we use the δτ x , δτ y , and δτ oi maps to infer the subsurface flows v.
We carry out the inversions of the travel-time maps with a noise level corresponding to what we would expect for an average over 100 supergranules measured for 24 hours each. To generate the input travel-time maps we start from the noise-subtraction travel-time maps ( Figures. 3 and 4) and add noise (computed from the simulations carried out in the flow-free reference model) with an amplitude reduced by 1/ √ 100. We use kernels computed in the Born approximation using the method of Birch and Gizon (2007) with erratum Birch, Gizon, and Burston (2011), and account for the holography Green's functions using the approach of . The vector-valued kernel functions K relate the flows in the interior to the travel-time maps, In the above equation, we have separate maps and kernel functions for each combination of ridge filter and frequency range. Figure 5 compares the measured travel-time shifts with the travel-time shifts that we would expect from Equation (9)  The label at the bottom of each column denotes the radial order of the filter (e.g. n = 0 is the f -mode, n = 1 is the p 1 -mode, etc.) and the factor by which the mode data has been scaled (e.g. ×3 means the data has been multiplied by a factor of three).
of the sensitivity kernels, for example the neglect of the sponge layers in the simulation. We discuss this issue in more detail in Sections §5 and 6. In order to discretize Equation (9) . Noise subtracted center-annulus travel-time differences δτ oi calculated using ridge filters and 0.25 mHz wide band-pass frequency filters for the supergranulation simulation. The label at the bottom of each column denotes the radial order of the filter (e.g. n = 0 is the f -mode, n = 1 is the p 1 -mode, etc.) and the factor by which the mode data has been scaled (e.g. ×3 means the data has been multiplied by a factor of three). The signal-to-noise ratio of the travel-time differences are generally greater than those of δτx (Figure 3).
where the sum over j is over a total of N basis functions φ j (z) and the sum over k is a two-dimensional inverse FFT on the grid that the travel-time maps are measured on (the simulation domain is periodic; as a result, the travel-time maps are periodic as well). Here we choose the basis functions φ j (z) to be one when z j < z ≤ z j+1 and zero otherwise. The points z i consist of a grid of 56 points that are equally spaced in acoustic depth and cover the range from 0.5  [9]). Each panel corresponds to a particular combination of frequency filter, radial order, and measurement geometry. The three images in each panel contain the forward model associated with the true flow (left), the measurement from the simulation (middle), and the residual (right), respectively. As discussed in the text, the forward model becomes less accurate at high frequencies. The maps have been smoothed with a gaussian filter (σ = ten pixels) to more clearly elucidate the differences, and the panes are scaled independently. Mm above the photosphere to 10 Mm below the photosphere. The grid spacing in k-space is h k = 0.063 rad Mm −1 .
We used the MCD approach (Jacobsen et al., 1999) to split the full inversion problem into small one-dimensional (z-only) inversion problems at each horizontal wavenumber. These problems we solved using RLS (e.g. Zhao, Kosovichev, and Duvall, 2001;Kosovichev, 1996, in the context of local helioseismology) with a regularization term given by the depth integral (at each horizontal wavevector k) of the quantity v 2 x + v 2 y + 10v 2 z ; the factor ten was chosen to reflect the preconception that vertical velocities are on average smaller than horizontal flows. Following Couvidat et al. (2005), we used a k-dependent regularization parameter with λ 2 = λ 2 1 + λ 2 2 k 2 , where (λ 1 , λ 2 ) are fixed parameters. As the regularization parameters are increased, the amplitude of the large-scale flow in the inversion result is reduced. We chose the parameters by inspection for each inversion based on a compromise between the desire to average out small-scale noise and the goal of retaining the large-scale flow features. Figure 6 shows a comparison of different regularization parameters applied to an inversion result for a portion of the v x velocity field.
The inversion results presented herein reflect an estimate of the mean flow averaged over 100 supergranules for a 24-hour window. Application to real data necessitates consideration of the effective averaging window relative to the time scale over which the coherent structure evolves; a typical supergranule lifetime is estimated to be a day or more (Worden and Simon, 1976;Hirzberger et al., 2008).

Results
Figures 7 and 8 show inversion results using all available n = 0, n = 1, and n = 2 measurements for (λ 2 1 , λ 2 2 /h 2 x ) = (3.24×10 −8 , 1.69×10 −4 ) Mm −1 (m s −1 ) −2 compared to the true v x and v z velocity components in the supergranule simulation, respectively. The structure of the v x -inverted fields in Figure 7 shares qualitative agreement with the true flow, although misses the signal at the top and bottom (below about z = −4 Mm) of the physical domain. The underestimation of flow magnitude in the lower layers is consistent with prior work (e.g. Zhao, Kosovichev, and Duvall, 2001). The v z inversion in Figure 8 fails to qualitatively capture the true flow throughout the physical domain. As we will discuss later, this is due to cross-talk effects between components of the velocity. We tested the importance of including the n = 2 measurements by recomputing the inversions shown in Figures 7 and 8 using only the n = 0 and n = 1 measurements. There was no discernible visual difference in the inversion results for the horizontal and vertical velocity due to exclusion of the n = 2 measurements, which indicates that the RLS technique is relatively insensitive to the addition of higher radial order measurements.
We tested the effect of noise in the measurements by recalculating the flow inversions using only the n = 0 and n = 1 measurements shown in Figures 3  and 4, which excludes the noisiest measurements (i.e. the whited-out boxes in the grid in Figures 3 and 4) that were included in the inversion calculations shown in Figures 7 & 8. There was no discernible visual difference in the inversion results for the horizontal and vertical velocity due to exclusion of the noisiest measurements, which supports the conclusion that the RLS technique is relatively insensitive to changes in the number and quality of measurements used in the inversion calculations. The inversion fails to qualitatively capture the true vz flow throughout the physical domain. As we will discuss later, this is due to cross-talk effects between signals.
We calculated center-quadrant travel-time differences from Equation (9) using the kernels and the inferred flow from Figures 7 and 8 and compared these to the measurements from the simulation. Figure 9 shows these comparisons and the calculated residual (difference between the two) for every combination of frequency filter, radial order, and measurement geometry. There is generally good correspondence between the perturbations in the n = 0 and n = 1 measurements; however, some of the δτ x and δτ y measurements, as well as all of the n = 2 measurements, contain significantly more noise. In some of the higher-frequency measurements, there is noticeable structure in the residual. This may be due to inaccuracies in the kernels (see Figure 5). We tested the sensitivity of the inversion to the high-frequency measurements by repeating the inversion using the same regularization parameters, but with the measurements at 4.25 mHz and 4.5 mHz removed. The main impact of removing these measurements was a reduction in the amplitude of the inferred vertical flow by about a factor of four.
Our findings are generally consistent with those of Zhao et al. (2007), who used the time-distance method and kernels computed from the ray approximation to perform inversions of simulated supergranule-scale convection. They were able to obtain reasonable inversion results only down to about 4 Mm below the surface, although they credit this to a limit in the largest annulus radius used. Here, the sensitivity functions for many of the travel-time measurements extend well below z = −4 Mm. This suggests that, with a sufficiently small noise level, it will be possible to detect flows below this depth (although with a spatial resolution that decreases with increasing depth). The structure of the horizontal OI, n=2 X, n=0 X, n=1 X, n=2 Y, n=0 Y, n=1 Y, n=2 Figure 9.: Travel-time measurements calculated from the forward model compared to measurements from the simulation. Each panel corresponds to a particular combination of frequency filter, radial order, and measurement geometry. The three images in each panel contain the forward model associated with the inversion result (left), the measurement from the simulation (middle), and the residual (right), respectively. All of the n = 0 and n = 1 measurements generally reflect the structure of the model supergranule flow; the n = 2 measurements reflect a significant noise contribution. For all radial orders, the δτ x and δτ y measurements contain more noise than the δτ oi measurements. The maps have been smoothed with a gaussian filter (σ = 10 pixels) to more clearly elucidate the differences, and the panes are scaled independently. The left-most panel shows a slice through the vx inversion result for reference. The figure demonstrates that the lobed nature of the vx inversion is due to the form of the kernels, which exhibit localized minima and maxima as a function of z. Furthermore, the sensitivity of the n = 2 kernels is generally less than the n = 0 and n = 1 kernels, which means that use of higher radial-order measurements in the inversion calculations will have decreasing effectiveness on the result.
velocity inversion results can be better explained by observing the form of the sensitivity kernels. Figure 10 contains a vertical slice through the v x inversion (left-most panel) and horizontally integrated δτ x kernels grouped by radial order and frequency, and normalized by the rms noise level. The n = 0 kernels have a single band of sensitivity with a maximum located in between 0 and -1 Mm. The n = 1 kernels have two regions of sensitivity which explains why the v x inversion has two local maxima in the vertical direction. The n = 2 kernels have several regions of sensitivity; however, inclusion of the n = 2 measurements did not significantly effect the performance of the inversions because the signal-to-noise ratio of the n = 2 measurements is about a factor of two lower than that of the n = 0 and n = 1 measurements ( Figure 4). Figure 11 shows the averaging kernel that relates the inferred v x to the true x-component of the velocity. The RLS inversion procedure is able to produce reasonably localized x-component averaging kernels with minimal leakage outside the target region. Zhao et al. (2007) discuss the influence of cross talk on the inversions for the vertical velocity. Specifically, they describe a region just below the photosphere where the v z inversion is of the incorrect sign. We find a similar phenomenon in our results, illustrated by the band of negative vertical velocity near z = 0 Mm (Figure 8). Here we will use the notation v ij to denote the contribution of the j component of the (known) model velocity to the inversion result for the i component of the velocity; computing each of the v ij involves the horizontal convolution of the appropriate averaging kernel with a particular component of the known model velocity (see, e.g. Jackiewicz, Gizon, and Birch, 2008, for examples). Figure 12 shows the contributions of each of the components of the model flow to the inferred vertical flow (i.e. v zx , v zy , v zz ). The contribution of the x-and y-components of the true flow to the inferred vertical flow are as large as the contribution from the true vertical flow. The sum of these velocities (v tot = v zx + v zy + v zz ) is qualitatively similar to the vertical flow inferred from the inversion. If the kernels were exactly correct and there was no noise then v tot would be identical to the inferred vertical flow. This computation highlights the importance of controlling the cross-talk in inversions for velocity (this has been achieved using OLA inversions by Jackiewicz, Gizon, and Birch, 2008;Švanda et al., 2011).

Discussion
We computed RLS inversions from travel-time measurements of simulated data created by numerically propagating waves through a simple supergranule-like flow field. The benefit of this approach is the opportunity for comparison with a known flow field for direct evaluation of the performance of the measurement and inversion techniques. Furthermore, the use of simulated data allows for isolation of complicating effects. For instance, here travel-time differences are due solely to the supergranule flow; travel-time differences calculated from observational solar data contain the effects of many more physical variables. We found that the inversion of the v x -and v y -velocity shared general features of the true flow reasonably well down to a depth of about 4 Mm below the surface. The inversion of the vertical velocity performed poorly throughout the domain, particularly near the surface where the computations produced a result that was of incorrect sign.
Underestimation of large-scale signal strength may be attributed to choice of regularization, which averages out noise but also causes some loss of flow ampli- tude. We tested the inversions over a broad range of regularization parameters, and the results presented here reflect the best scenarios.
We found the inversions, and hence the averaging kernels, to be relatively irresponsive to an increase in the number of input travel-time maps or to exclusion of noisy travel-time maps. This is a credit to the capability of the RLS technique in minimizing the effects of noise; however, the significant error in the inversion results suggests that the forward modeling is inadequate.
The sensitivity kernels relate variations in travel-time shifts to variations in the solar model; capturing this sensitivity accurately is dependent on making appropriate assumptions regarding the physical state (in our case, the simulated state) and the measurement procedure (e.g. Gizon and Birch, 2002;Birch, Kosovichev, and Duvall, 2004). As a test of the kernels, we convolved the resulting averaging kernels with the components of the true supergranule flow field, and compared these results with the respective inversion results. We found that the convolutions contained qualitatively the same structure (and hence comparable errors) as the inversion results.
There is very little sensitivity to flows above the photosphere because only waves of frequency comparable to the acoustic cutoff frequency (about 5.5 mHz) sample this region. The calculation of the kernels uses a zero Lagrangian pressure perturbation upper boundary condition, an assumption which is likely not accurate for high-frequency waves. Furthermore, the simulation uses absorbing sponges at the top and bottom boundaries, which may affect the wave dynamics significantly differently from what is accounted for in the kernel calculations. Perhaps for these reasons, the inversions are unable to reproduce the flow structure near the top of the computational domain. Other constraints on the forward modeling and approximations used in the design of the sensitivity kernels may need to be considered in order to improve the performance of flow inversions.
In conclusion, our findings indicate that errors in the RLS inversions result from the combined effects of cross talk between signals, choice of regularization parameters, design of sensitivity kernels, and inconsistencies among boundary conditions, all of which are open areas of research requiring further study. We note that improvement in performance, particularly in regard to cross-talk effects, has been seen in the optimally localized average (OLA) approach to reduce the size of the side lobes on the averaging kernels and limit spatial leakage (Jackiewicz, Gizon, and Birch, 2008;Švanda et al., 2011;Jackiewicz et al., 2012).