Supergranulation and multiscale flows in the solar photosphere
Key Words.:
Sun: photosphere – Sun: interior – Convection – Turbulence – InstabilitiesThe Sun provides us with the only spatially wellresolved astrophysical example of turbulent thermal convection. While various aspects of solar photospheric turbulence, such as granulation (oneMegameter horizontal scale), are well understood, the questions of the physical origin and dynamical organization of largerscale flows, such as the 30Megameters supergranulation and flows deep in the solar convection zone, remain largely open in spite of their importance for solar dynamics and magnetism. Here, we present a new critical global observational characterization of multiscale photospheric flows and subsequently formulate an anisotropic extension of the BolgianoObukhov theory of hydrodynamic stratified turbulence that may explain several of their distinctive dynamical properties. Our combined analysis suggests that photospheric flows in the horizontal range of scales between supergranulation and granulation have a typical vertical correlation scale of 2.5 to 4 Megameters and operate in a strongly anisotropic, selfsimilar, nonlinear, buoyant dynamical regime. While the theory remains speculative at this stage, it lends itself to quantitative comparisons with future highresolution acoustic tomography of subsurface layers and advanced numerical models. Such a validation exercise may also lead to new insights into the asymptotic dynamical regimes in which other, unresolved turbulent anisotropic astrophysical fluid systems supporting waves or instabilities operate.
1 Introduction
The solar photosphere is the stage of many spectacular (magneto)hydrodynamic phenomena and it provides us with a unique observationally wellresolved example of strongly nonlinear thermal convection, one of the most common fluid instabilities and transport processes encountered in nature and astrophysics (Kupka et al. 2007). While some aspects of solar thermal turbulence, such as granulation, are now well understood (Nordlund et al. 2009), we still lack definitive answers to many important questions, such as how the turbulence organizes on large scales, and how it interacts with and amplifies magnetic fields or transports quantities such as angular momentum or magnetic flux (Miesch 2005; Charbonneau 2014).
The most direct observational characterization of flows connected to the solar convection zone at scales larger than the granulation has been achieved through measurements of Dopplerprojected velocities at the photospheric level. In particular, even a simple visual inspection of Doppler images of the quiet Sun clearly reveals the pattern of supergranulation flow “cells”, whose trademark signature is a peak around (35 Megameters, Mm) in the sphericalharmonics energy spectrum of Dopplerprojected velocities (Hathaway et al. 2000; Williams et al. 2014; Hathaway et al. 2015). The physical origin of this supergranulation is widely debated (Rieutord & Rincon 2010). In particular, while the idea that supergranulationscale flows may just be a manifestation of some form of thermal convection has a long history, it has often been dismissed due to the seeming lack of photometric intensity contrast at the same scales (Langfellner et al. 2016). Besides, there is as yet no general consensus on local helioseismic estimates of the amplitude, depth and structure of subsurface flows on this scale (Nordlund et al. 2009; Rieutord & Rincon 2010; Gizon et al. 2010; Švanda et al. 2013; DeGrave et al. 2014; Švanda 2015) – or in fact on any scale larger than that (see, e.g., Hanasoge et al. 2012; Gizon & Birch 2012; Hanasoge et al. 2016; Greer et al. 2015; Toomre & Thompson 2015; Greer et al. 2016). Numerical simulations of the problem have, until recently, also been quite limited. While their results have not led to clearcut results and conclusions (see reviews by Miesch 2005; Nordlund et al. 2009; Rieutord & Rincon 2010), many of them suggest that meso to supergranulationscale flows have a convective (buoyant) origin (e.g., Rincon et al. 2005; Lord et al. 2014; Cossette & Rast 2016).
In this article, we attempt to make further progress on this problem using the highquality observations of the solar photosphere provided by the Solar Dynamics Observatory (SDO) and a new theoretical analysis. In Sect. 2, we present a global observational analysis of multiscale photospheric vector flows reconstructed from quasifull disc Doppler and photometric data from SDO using a tracking technique. This analysis subsequently leads us in Sect. 3 to formulate a theory of anisotropic turbulent convection that may explain several distinctive dynamical properties of these flows. As further argued in Sect. 4, the combined results consistently suggest that photospheric flows in the horizontal range of scales in between supergranulation and granulation operate in a strongly anisotropic, selfsimilar, buoyant dynamical regime. Connections between these results and earlier work, as well as future perspectives, are discussed in Sect. 5. The main text of the article is focused on results. The technical details of the data processing and analysis are provided in the two Appendices.
2 Observations
Our observational analysis is based on 24 hours of uninterrupted highresolution whitelight intensity and Doppler observations of the entire solar disc by the HMI instrument aboard the SDO satellite (Scherrer et al. 2012; Schou et al. 2012). The data was obtained on October 8, 2010 starting at 14:00:00 UTC, one of the quietest periods of solar activity since the launch of SDO. Our analysis is distinct and independent of the recent studies by Williams et al. (2014); Langfellner et al. (2015); Hathaway et al. (2015), and both confirms and extends some of their results.
2.1 Data analysis
2.1.1 Velocity field reconstruction
Our first objective is to reconstruct the three components of the photospheric vector velocity field in the angular degree range from available observations. To this end, we perform a Coherent Structure Tracking analysis (CST, Roudier et al. 2012) of photometric structures such as granules. This technique provides us with the projection of the photospheric velocity field onto the plane of the sky/CCD matrix (Fig. 1). We can then combine the results with Doppler data, which provides the outofplane component of the velocity field, to calculate the full vector velocity field at the surface in solar spherical coordinates,
(1)  
(2)  
(3)  
where is the heliographic latitude of the central point of the solar disc at the time of observation. The reduction and filtering of the raw photometric and Doppler data required to build a consistent data set for is described in detail in App. A.
2.1.2 Sphericalharmonics analysis
The reconstructed vector field is expanded in terms of vector spherical harmonics as
(4) 
where , , and . The spectral coefficients describe the radial flow component, and the coefficients and describe the spheroidal (divergent) and toroidal (vortical) parts of the horizontal flow components, respectively. Global maps of the horizontal divergence and radial vorticity at the surface are obtained by mapping back and to physical space. The corresponding onedimensional radial, spheroidal and toroidal energy spectra are given by
(5)  
(6)  
(7) 
respectively. The prefactors in and result from the definition of the spheroidal and toroidal vector spherical harmonics basis vectors.
The technical aspects of the harmonictransforms procedure of SDO/HMI data are described in detail in App. B. Here, it suffices to mention that the data must first be apodized because we only see one side of the Sun, and we must eliminate nearlimb regions where the CST analysis is not reliable. Besides, as will be shown below, the radial component of the flow is very weak and can only be tentatively determined using data close to the disc center to limit contamination by the intrinsic algorithmic noise of the CST in the deprojection process. We therefore use an apodizing window with a opening heliocentric angle from the disc center to compute the spectra of the radial and lineofsight (Doppler) velocity fields, and a window with a opening angle to compute the spectra of the horizontal velocity field. The general mathematical definition of these windows is provided in Eq. (37).
2.2 Results
2.2.1 Velocity field maps
Figure 2 shows maps of the horizontal divergence field and vertical vorticity field derived from a single snapshot of the 30minuteaveraged horizontal surface velocity field using the sphericalharmonics spheroidaltoroidal decomposition into divergent and vortical motions. The pattern is dominated by supergranulation, whose divergent morphology is revealed in the closeup panel using a lineintegral convolution technique highlighting horizontalflow Lagrangian trajectories (Cabral & Leedom 1993). The vortical component of the flow is significantly weaker than the divergent component, as indicated by the associated rate of strains reported in the color bars (and further demonstrated in the spectral analysis presented below). Figure 3 shows a local map in the same disccenter zone of the 24haveraged , superimposed with horizontal flows averaged over the same time scale. There is a noticeable correlation between the weak upflows (downflows) and horizontal divergences (convergences) at the supergranulation scale.
2.2.2 Energy spectra
Figure 4 shows the sphericalharmonics radial, spheroidal and toroidal energy spectra , , and of the 30minuteaveraged velocityfield components. The horizontal divergent component is energetically dominant at all scales probed: the corresponding spheroidal spectrum exhibits a clear peak at the supergranulation scale , corresponding to a typical r.m.s. flow velocity of . The r.m.s. radial flow is at . The spectra of the different flow components behave differently for : the spheroidal and radial spectra respectively decay and increase with decreasing scale, while the toroidal spectrum is almost flat. The exact shapes of the weaker and noisier and are moderately sensitive to the processing and averaging of the data, but the general trend described here is robust (as explained in App. A.2, it is very difficult to separate the Doppler contributions of the largescale spectral tail of granulation from those of the radial component of the largescale motions otherwise detected by the CST. Consequently, and while we did our best to remove the largescale contribution of granulation from the Doppler data, our results for in the range still probably slightly overestimate the actual amplitude of the radial component of flows not associated with granulation. We also find that lies somewhere between and power laws, depending on the exact filtering procedure applied to the Doppler signal). The spectral falloff at is due to the intrinsic 2.5–5 Mm CST spatial resolution down to which the vectorfield reconstruction was performed. The different power laws that these apparently selfsimilar spectra follow for will be discussed in Sects. 3.3 and 4.3.
2.2.3 Vertical correlation scale
Finally, we estimate the typical vertical correlation scale of these subsonic flows as a function of horizontal scale (or of the angular degree ). This can be obtained from the mass conservation relation
(8) 
where and are the horizontal and vertical velocity fluctuations (because the medium is stratified, here serves as an a priori proxy notation for either the radial wavelength of the perturbations or the local density scale height, whichever one is smaller). We therefore define
(9) 
and plot this quantity in Fig. 4 (inset) for . Remarkably, the typical vertical correlation scale does not vary significantly in the range ( Mm), that is Mm. Considering the uncertainties in the determination of discussed in Sect. 2.2.2, this scale estimate should probably be considered as an upper bound rather than as an exact value. Overall, this scale appears to be comparable with the density scale height below the surface, but is somewhat deeper than the granulation thermal boundary layer (Nordlund et al. 2009).
2.3 Main conclusions
The observational analysis presented above leads us to the following conclusions.

The flow in the range is characterized by divergent (convergent) horizontal flows correlated with upflows (downflows). This predominantly celllike morphology of the flow suggests that the dynamics in that range of scales is dominated by buoyancy forces.

The ratio of to suggests that motions are strongly anisotropic and strongly feel the effects of stratification, with the vertical correlation scale of order the density scale height.

The significant energy content of horizontal motions and their spectral break at the supergranulation scale, followed by the apparent powerlaw decay of their spectrum at smaller horizontal scales, suggest that supergranulation corresponds to the largest buoyancydriven scale of the system, below which some form of nonlinear cascade takes place.
3 Theoretical considerations
Let us now attempt a theory of anisotropic turbulent convection that could explain these observations. As argued by Rincon (2007), a good starting point is the BolgianoObukhov (BO59, see Oboukhov 1959; Bolgiano 1962) phenomenology, which is based on the following assumptions: (i) a constant spectral flux of buoyancy fluctuations follows from the temperature/energy equation, (ii) a balance between the inertial and buoyancy terms in the momentum equation, and (iii) isotropy of motions. BO59 predicts a thermal spectrum, and a steep velocity spectrum. For the Sun, the latter is ruled out by Fig. 4 for . This is not very surprising because, while the first two BO59 assumptions make sense in the present context, the isotropy assumption must manifestly be false because (see Sect. 2.2.3).
More generally, we point out that the isotropic BO59 scaling laws are not observed either in simulations of RayleighBénard convection with orderunity vertical to horizontal aspect ratio (Rincon 2006, 2007; Kumar et al. 2014). While the conclusions of these previous studies apply to the regime (where is an isotropic wave number), they do not appear to be directly applicable to the strongly anisotropic regime characteristic of the astrophysical problem investigated in the present paper. In any case, in what follows, we only use the phenomenology of the original isotropic BO59 theory as a physical motivation to derive distinct anisotropic scaling laws for the largely unexplored “long horizontal wavelength” regime .
3.1 Theoretical framework
An anisotropic generalization of BO59 may be derived from the following set of dynamical equations valid in the Boussinesq and anelastic limits:
(10)  
(11)  
(12) 
Here is the mean density, denotes velocity perturbations, denotes pressure fluctuations divided by , is the actual temperature in the Boussinesq limit, and the potential temperature in the anelastic limit, is its horizontally averaged part and its fluctuations around this average, is the kinematic viscosity, is the thermal diffusivity ( in the Sun), and is the acceleration of gravity.
3.2 Phenomenology of anisotropic turbulent convection
The solar convection zone underneath the thin thermal boundary layer at the surface is in the strongly nonlinear convection regime, so its thermodynamic profile is very close to adiabatic and is almost independent of . We may then argue that the thermal injection term on the r.h.s. of Eq. (11) is small compared to the nonlinear advection term on the l.h.s.. Similarly to the Kolmogorov phenomenology in hydrodynamics (K41, see, e.g., Davidson 2004), which postulates a constant flux of kinetic energy in the inertial range (, the dynamics in this “thermalinertial” regime is characterized by a constant flux of variance from scale to scale:
(13) 
where is the (time) average dissipation rate of thermal fluctuations (equal to their injection rate) and denotes increments of fluctuations between two points separated by a horizontal separation vector at a constant radial coordinate. Equation (13) effectively corresponds to BO59 assumption (i). This can, in fact, also be formalized on the basis of Yaglom’s exact law for dissipative scalars (Monin & Yaglom 1971; Rincon 2006), which Eq. (11) obviously satisfies if . Similarly to BO59 assumption (ii), we seek a buoyancydominated dynamical regime balancing the inertial and buoyancy terms in Eq. (10), but we drop the isotropy assumption, expecting instead that . The result of this is
(14) 
The small prefactor on the r.h.s. accounts for the partial cancellation of buoyancy by the pressure gradient in the vertical direction, as a result of the horizontal redistribution of buoyancy work by pressure forces (this prefactor also appears in the dispersion relation of gravity waves in the limit ).
3.3 Scaling laws
Equations (13)(14), combined with the mass conservation relation Eq. (8), lead to the following selfsimilar scaling relations,
(15)  
(16)  
(17) 
As a consistency check, we can see that the isotropic BO59 prediction , is recovered if we set . Here however, we are instead interested in a strongly stratified regime in which the typical vertical scale height of the system is much smaller than . We therefore treat as a simple parameter equal to in Eqs. (15)(17), as also suggested by our observational result of Sect. 2.2.3. Doing this results in the following powerlaw energy spectra as functions of horizontal wave number ,
(18)  
(19)  
(20) 
We note that in the regime , relevant to the astrophysics problem investigated here, these scalings are distinct from the classical isotropic BO59 prediction.
4 Comparison of theory with observations
4.1 Selfsimilar buoyant dynamics
Physically, the theory proposed above describes a situation in which buoyant thermal fluctuations stochastically injected through the surface boundary layer drive nonlinear anisotropic convective motions (plumes) in the adiabatic layer. By mixing these thermal fluctuations, turbulent velocity fluctuations drive a horizontal cascade of buoyancy down to thermal dissipative scales, resulting in statistically selfsimilar buoyant dynamics. This picture appears to be very much in line with the observed phenomenology of flows at scales larger than the granulation, notably the hierarchical, bubbling recurrent dynamics of “trees and families of fragmenting granules” documented by Roudier et al. (2003, 2009, 2016). The dynamics are ’quasi2D’ in the sense that there is no refinement of the vertical scale as the horizontal cascade proceeds, but it is 3D in the sense that finite is required to ensure mass conservation.
4.2 Spectral amplitudes vs. theoretical estimates
We now test for consistency between the observed spectral amplitudes and our theoretical scaling estimates derived by assuming that the observed fluctuations have a convective origin. A fit of the observed to Eq. (19) in the range for gives
(21) 
On the other hand, a theoretical dimensional estimate of the same quantity can be obtained directly in terms of the solar luminosity and internal structure of the Sun. The convective flux as a function of horizontal scale, as estimated from Eqs. (15)(17), is
(22) 
where is the specific heat at constant volume. We see that this flux increases with decreasing horizontal scale (the horizontal kinetic energy, on the other hand, increases with increasing horizontal scale in this regime. This is a consequence of the different anisotropic scalings for and , and suggests that supergranulationscale convection, while dominant in terms of horizontal kinetic energy, does not transport a significant fraction of the solar flux). However, this scaling theory is only applicable for . As becomes of the order of , there should be a crossover to either the standard isotropic Bolgiano regime , , or directly to the Kolmogorov regime , (Kumar et al. 2014). In both cases, the cascade proceeds in both horizontal and vertical directions, with for BO59 or for K41, that is, the convective flux decreases with decreasing scale . We therefore argue that the convective flux peaks at the isotropization scale , and require that at this scale be of the order of the total solar flux,
(23) 
where is the solar luminosity. Using , , (Stix 2004) and in Eq. (23), we obtain the theoretical estimate
(24)  
which, given the approximate nature of the calculation and the neglect of geometric prefactors in Eqs. (15)(17) and Eqs. (18)(20), is relatively close to our observational estimate in Eq. (21) (this cascade time scale also appears to be similar to the typical plume thermal diffusion time scale in the Sun (Miesch 2005; Cossette & Rast 2016), that is, the “turbulent” diffusivity at scale is of the order of the molecular diffusivity ). Overall, the observed amplitudes therefore seem consistent with a convective origin of the turbulence. We note, finally, that if this estimate holds, then Eq. (17) suggests relative temperature fluctuations of approximately at supergranulation scales in the first few Mm below the surface, consistent with helioseismic inferences (e.g., Duvall Jr et al. 1997).
4.3 Spectral power laws
The horizontal spectral scaling laws predicted by our theory for the vertical and horizontal velocities, (Eq. (18)) and (Eq. (19)), are broadly consistent with our observations, as shown in Fig. 4 by the power laws in the relevant range of scales. Different tentative scalings and , derived by naively postulating K41 scalings for the horizontal velocity field and using the mass conservation Eq. (8) with to obtain the scaling of the vertical velocity component, are shown for comparison (dotted lines). Such scalings cannot be ruled out given the closeness of their spectral exponents to those of our theory and the aforementioned sensitivity in the determination of to the details of data processing. However, we see no reasonable physical argument as to why or how such an “anisotropic K41” regime would actually be realized in this environment, in contrast with the nonlinear buoyancy regime described above which naturally produces divergent horizontal motions at all scales (which are manifest in Figs. 2 and 3). At this stage, we do not have a theory for the scalings of the energy spectrum of the weaker horizontal toroidal velocity field, which seems to be only weakly coupled to the spheroidal divergent dynamics.
Finally, at even larger scales (), beyond the supergranulation spectral break, our observational analysis shows a regime change to an intriguing . This might be one of the first observations of a largescale preinjection spectrum predicted by Saffman (1967). In theory, this spectrum is established as a result of longrange pressure correlations in turbulent flows in which fluid “eddies” have a random distribution of linear momentum (Davidson 2004), as opposed to a random distribution of angular momentum, which would instead result in a spectrum (Batchelor & Proudman 1956). In the solar photosphere, such eddies could perhaps be associated with supergranulation cells and smallerscale convective motions.
5 Discussion
We have presented a set of consistent arguments suggesting that turbulent flows at scales larger than granulation are a manifestation of statistically selfsimilar, strongly nonlinear convection in the bulk adiabatic layer below the solar surface. While the classic idea that supergranulation has a convective origin has usually been dismissed due to the seeming lack of photometric intensity contrast at the corresponding scales at the surface (Langfellner et al. 2016), helioseismology suggests that relative temperature fluctuations of approximately a few percent are present underneath the surface at such scales (Duvall Jr et al. 1997), in line with our calculation in Sect. 4.2. Besides, all existing largescale simulations, including those with radiative transfer that do not display any “meso” or “super” scale intensity contrast at the surface, exhibit a strong buoyancy driving and flows at such scales in the bulk of the convective layer (e. g., Rincon et al. 2005; Nordlund et al. 2009; Bushby & Favier 2014; Lord et al. 2014; Cossette & Rast 2016, see also Rieutord & Rincon (2010), Fig. 12). Importantly, in these simulations, these scales are significantly larger than those of the most unstable linear modes of the system at rest, and the statistical order at such scales emerges dynamically in the nonlinear regime. If there are indeed significant thermal fluctuations associated with supergranulationscale convection in the first few Mm below the surface, they may be blanketed by the thinner thermal granulation boundary layer, or could be optically thick due to the extreme dependence of opacities on temperature (Nordlund et al. 2009). Our theory also suggests that anisotropic supergranulationscale convection, although it generates strong horizontal flows, does not significantly contribute to the transport of thermal energy due to the weakness of the radial component of the velocity field at such scales.
We do not have a quantitative answer as to what determines the spectral supergranulation break, and can only offer directions. Our theory predicts that thermal fluctuations increase with scale as in the nonlinear convection regime. However, there is a physical limit to how large such fluctuations can be, which could determine the maximum buoyant scale at which the selfsimilar scaling should break down. This limit should be related in some way to the maximum entropy fluctuations that can be injected into the adiabatic layer and, therefore, to the details of granulation and of the superadiabatic thermal boundary layer at the surface. This conclusion appears to be in line with the recent numerical simulations of this problem by Cossette & Rast (2016), which show that the thickness of the boundary layer has a strong effect on the energycontaining scale of the turbulent energy spectra. However, other physical explanations (e.g., Featherstone & Hindman 2016) cannot be ruled out at this stage.
Based on our experience with this data, the observational analysis and comparison with theory presented above are close to the limit of what is achievable with a combination of surface tracking and direct Doppler measurements given their disparate natures. Future progress will probably come from acoustic tomography (Toomre & Thompson 2015) and advanced numerical models. While somewhat speculative at this stage, our theory, including Eqs. (18)(20), may soon be testable using such tools. An important question in this respect asks to what extent pristine powerlaw scalings derived from simple dynamical arguments can be realized in systems such as the photospheric transition region, where a variety of physical processes become intertwined.
The qualitative implications of such analyses could be wider than the solar context. They may notably provide us with fundamental insight into the structure of anisotropic turbulence in general (e.g., Nazarenko & Schekochihin 2011), as well as into turbulence in more distant, unresolved astrophysical systems supporting anisotropic waves or instabilities, such as stellar interiors (Zahn 2008; Miesch & Toomre 2009), galaxy clusters (Zhuravleva et al. 2014) and accretion discs (Walker et al. 2016).
References
 Batchelor & Proudman (1956) Batchelor, G. K. & Proudman, I. 1956, Phil. Trans. R. Soc. Lond., Ser. A, 248, 369
 Bolgiano (1962) Bolgiano, R. 1962, J. Geophys. Res., 67, 3015
 Bushby & Favier (2014) Bushby, P. J. & Favier, B. 2014, A&A, 562, A72
 Cabral & Leedom (1993) Cabral, B. & Leedom, L. C. 1993, in Proceedings of the 20th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH ’93 (New York, NY, USA: ACM), 263
 Charbonneau (2014) Charbonneau, P. 2014, Annu. Rev. Astron. Astrophys., 52, 251
 Cossette & Rast (2016) Cossette, J.F. & Rast, M. P. 2016, [arXiv:1606.04041]
 Davidson (2004) Davidson, P. A. 2004, Turbulence (Oxford University Press)
 DeGrave et al. (2014) DeGrave, K., Jackiewicz, J., & Rempel, M. 2014, ApJ, 788, 127
 Duvall Jr et al. (1997) Duvall Jr, T. L., Kosovichev, A. G., Scherrer, P. H., et al. 1997, Solar Phys., 170, 63
 Featherstone & Hindman (2016) Featherstone, N. A. & Hindman, B. W. 2016, ApJ in press [arXiv:1609.05153]
 Gizon & Birch (2012) Gizon, L. & Birch, A. C. 2012, Proc. Natl. Acad. Sci., 109, 11896
 Gizon et al. (2010) Gizon, L., Birch, A. C., & Spruit, H. C. 2010, Annu. Rev. Astron. Astrophys., 48, 289
 Greer et al. (2015) Greer, B. J., Hindman, B. W., Featherstone, N. A., & Toomre, J. 2015, ApJ, 803, L17
 Greer et al. (2016) Greer, B. J., Hindman, B. W., & Toomre, J. 2016, ApJ, 824, 128
 Hanasoge et al. (2016) Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Annu. Rev. Fluid Mech., 48, 191
 Hanasoge et al. (2012) Hanasoge, S. M., Duvall, T. L., & Sreenivasan, K. R. 2012, Proc. Natl. Acad. Sci., 109, 11928
 Hathaway (1992) Hathaway, D. H. 1992, Solar Phys., 137, 15
 Hathaway et al. (2000) Hathaway, D. H., Beck, J. G., Bogart, R. S., et al. 2000, Solar Phys., 193, 299
 Hathaway et al. (2015) Hathaway, D. H., Teil, T., Norton, A. A., & Kitiashvili, I. 2015, ApJ in press, [arXiv:1504.01422]
 Howe & Thompson (1998) Howe, R. & Thompson, M. J. 1998, A&AS, 131, 539
 Kumar et al. (2014) Kumar, A., Chatterjee, A. G., & Verma, M. K. 2014, Phys. Rev. E, 90, 023016
 Kupka et al. (2007) Kupka, F., Roxburgh, I., & Chan, K. L., eds. 2007, IAU Symposium, Vol. 239, Convection in Astrophysics (IAU S239)
 Langfellner et al. (2016) Langfellner, J., Birch, A. C., & Gizon, L. 2016, A&A, 596, A66
 Langfellner et al. (2015) Langfellner, J., Gizon, L., & Birch, A. C. 2015, A&A, 581, A67
 Lord et al. (2014) Lord, J. W., Cameron, R. H., Rast, M. P., Rempel, M., & Roudier, T. 2014, ApJ, 793, 24
 Miesch (2005) Miesch, M. S. 2005, Living Rev. Solar Phys., 2, 1
 Miesch & Toomre (2009) Miesch, M. S. & Toomre, J. 2009, Annu. Rev. Fluid Mech., 41, 317
 Monin & Yaglom (1971) Monin, A. S. & Yaglom, A. M. 1971, Statistical Fluid Mechanics: Mechanics of Turbulence (MIT Press)
 Nazarenko & Schekochihin (2011) Nazarenko, S. V. & Schekochihin, A. A. 2011, J. Fluid Mech., 677, 134
 Nordlund et al. (2009) Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Living Rev. Solar Phys., 6, 2
 Oboukhov (1959) Oboukhov, A. M. 1959, Dokl. Akad. Nauk. SSR, 125
 Rieutord & Rincon (2010) Rieutord, M. & Rincon, F. 2010, Living Rev. Solar Phys., 7, 2
 Rieutord et al. (2007) Rieutord, M., Roudier, T., Roques, S., & Ducottet, C. 2007, A&A, 471, 687
 Rincon (2006) Rincon, F. 2006, J. Fluid Mech., 563, 43
 Rincon (2007) Rincon, F. 2007, in IAU Symposium, Vol. 239, Convection in Astrophysics, ed. F. Kupka, I. Roxburgh, & K. L. Chan, 58–63
 Rincon et al. (2005) Rincon, F., Lignières, F., & Rieutord, M. 2005, A&A, 430, L57
 Roudier et al. (2003) Roudier, T., Lignières, F., Rieutord, M., Brandt, P. N., & Malherbe, J. M. 2003, A&A, 409, 299
 Roudier et al. (2016) Roudier, T., Malherbe, J. M., Rieutord, M., & Frank, Z. 2016, A&A, 590, A121
 Roudier et al. (2009) Roudier, T., Rieutord, M., Brito, D., et al. 2009, A&A, 495, 945
 Roudier et al. (2012) Roudier, T., Rieutord, M., Malherbe, J. M., et al. 2012, A&A, 540, A88
 Roudier et al. (2013) Roudier, T., Rieutord, M., Prat, V., et al. 2013, A&A, 552, A113
 Saffman (1967) Saffman, P. G. 1967, J. Fluid Mech., 27, 581
 Schaeffer (2013) Schaeffer, N. 2013, Geochemistry, Geophysics, Geosystems, 14, 751
 Scherrer et al. (2012) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Solar Phys., 275, 207
 Schou & Brown (1994) Schou, J. & Brown, T. M. 1994, A&AS, 107, 541
 Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Solar Phys., 275, 229
 Stix (2004) Stix, M. 2004, The Sun : an introduction (Springer)
 Toomre & Thompson (2015) Toomre, J. & Thompson, M. J. 2015, Space Science Reviews, 196, 1
 Švanda (2015) Švanda, M. 2015, A&A, 575, A122
 Švanda et al. (2013) Švanda, M., Roudier, T., Rieutord, M., Burston, R., & Gizon, L. 2013, ApJ, 771, 32
 Walker et al. (2016) Walker, J., Lesur, G., & Boldyrev, S. 2016, MNRAS, 457, L39
 Williams et al. (2014) Williams, P. E., Pesnell, W. D., Beck, J. G., & Lee, S. 2014, Solar Phys., 289, 11
 Zahn (2008) Zahn, J.P. 2008, in IAU Symposium, Vol. 252, The Art of Modeling Stars in the 21st Century, ed. L. Deng & K. L. Chan, 47
 Zhuravleva et al. (2014) Zhuravleva, I., Churazov, E., Schekochihin, A. A., et al. 2014, Nature, 515, 85
Acknowledgements.
This work is dedicated to the memory of JeanPaul Zahn, a pioneer in astrophysical fluid dynamics research whose passion and work have been an invaluable source of inspiration to us. We thank Nathanaël Schæffer for his assistance with the SHTns library, the SDO/HMI data provider JSOC, the SDO/HMI team members, in particular P. Scherrer, S. Couvidat and J. Schou for sharing information on the calibration and removal of systematics. This work was granted access to the HPC resources of CALMIP under allocation 2011[P1115]. We also thank N. Renon for his assistance with the CST code parallelization. The work of AAS was supported in part by grants from the UK STFC and EPSRC. FR and AAS thank the Wolfgang Pauli Institute, Vienna, for its hospitality.Appendix A Photospheric velocity field reconstruction
a.1 Data reduction
In order to obtain needed for the reconstruction of the three components of the photospheric velocity field in Eqs. (1)(3), we used raw HMI photometric and Doppler images with a pixel size () every .
a.1.1 Reorientation and derotation of images
All the images were first reoriented using the SDO CROTA2 fits header keyword to have the solar north pole pointing ’up’ ( corresponds to the northern hemisphere), and derotated’ and aligned to avoid rotational smearing. We used the following rotation profile adjusted from the raw Doppler data averaged over the 24 hours of observations:
(25) 
where is the latitude.
a.1.2 Coherent Structure Tracking (CST) analysis
Starting with a sequence of intensity images, the CST algorithm (Rieutord et al. 2007; Roudier et al. 2012, 2013) determines the projection (,) in the CCD plane of the photospheric Eulerian velocity field by following the trajectories of many individual intensity structures such as granules during their entire “life” (the period of time, typically a few minutes, during which a structure remains a coherent, single object that does not split or merge). The average velocity of each structure is computed from its total displacement over this time interval, and supersonic velocities (larger than 7 ) are rejected. The velocities obtained this way are irregularly spaced on the field of view, and the last part of the procedure consists in approximating the results by the best differentiable Eulerian field fitting the data, using a multiresolution wavelet analysis (Rieutord et al. 2007). Applying this procedure to the SDO/HMI data set, we obtained a sequence of velocityfield maps with a temporal resolution of 30 min and a spatial resolution of 2.5 Mm (3.5”, or 7 HMI pixels, corresponding to a maximum angular degree ). We further removed the velocity signal associated with the satellite motion
(26)  
(27) 
where is the distance from the satellite to the Sun, and and are the westwards and northwards components of the satellite peculiar velocity (given in the data file’s headers). We finally found it necessary to filter out all the signal at spherical harmonics , as the latter is polluted by systematic largescale velocityfield residuals associated with imperfect corrections of the rotation, satellite motions, plages regions, and artificial largescale drifts imprinted by the CST analysis itself far from disc center (Roudier et al. 2013). This filtering was performed in the sphericalharmonics space by means of the harmonictransform machinery described App. B.
a.1.3 Doppler data analysis
Dopplergrams map the outofplane (lineofsight, l.o.s.) component of the photospheric velocity field. The SDO/HMI convention is that is taken as positive when the flow is away from the observer, so that the outofplane velocity towards the observer is (Fig. 1). This equality is only approximate because of the inclination between the actual l.o.s and the unit vector due to the finite distance between the Earth and the Sun. This effect is taken into account to eliminate the global rotational satellite motion signals, but is otherwise negligible. We first substracted the rotational Doppler shift
(28) 
where is given in Eq. (25), and the Doppler shift associated with the satellite motion
(29) 
where . The images were subsequently ’derotated’ in the same way as the whitelight images and further corrected from a polynomial radial limbshift function adjusted from ring averages of two hours of data taken at different heliocentric angles from disc center ,
(30) 
where and
(31) 
Since the velocityfield maps derived from the CST have a 2.5 Mm spatial resolution, we then downsampled the the Dopplergrams to , keeping only one point out of seven in each direction, and averaged over 30min periods to obtain an effective sequence of maps of , with the same spatial and temporal resolution as that of the maps derived from CST.
a.2 Noise analysis and  filtering
While the raw CST signal is not sensitive to solar oscillations or granulation, it is contaminated by an undesirable intrinsic algorithmic noise that can, to some extent, be removed via  filtering using the harmonictransform machinery described in App. B. The  diagram of the component of the raw CST velocity field computed from 24 hours of data is shown in Fig. 5 (top – the component behaves similarly). The main signal (red bump) is contaminated by a whitenoiseintime component in the form of a vertical white stripe. To remove this noise, we filtered out all the data at frequencies larger than those given by the filtering envelope (also shown in the figure):
(32) 
The  diagram of the Doppler signal (first downgraded to the CST resolution) for the same range of frequencies is shown in Fig. 5 (bottom ; the highfrequency part, not shown here, contains the usual mode ridges). The spectral support of the turbulent Doppler signal appears to be somewhat different from that of the CST. In particular, there is a fairly large, relatively highfrequency component in the range (in white in the figure), which we attribute to the largescale spectral tail of granulation and which is not captured by the CST. Note that the decrease of this signal at the smallest scales represented is due to the coarsegraining associated with the downgrade of the Doppler maps to the CST resolution.
The presence of this extra signal component in the Doppler data but not in the CST implies that the two raw signals cannot naively be mixed in the 3Dvelocityfield reconstruction in Eqs. (1)(3), as this would lead to spurious crosstalk between the different field components. In our view, the most consistent way to proceed is, therefore, to apply the same  filtering envelope (32) to both Doppler and CST fields. While it is imperfect and does not totally remove the granulation tail contribution from the largescale signal, this procedure focuses the analysis on largescale motions inferred from the CST, limits the risk of crossanalysing signals of different physical origin, and makes it possible to extract the weak supergranulationscale signal from the Doppler data out of the strong largescale tail of the surface granulation dynamics. This filtering also entirely removes the solar oscillation signal from the Doppler data.
Finally, we found it necessary to discard the residual signal in the spherical harmonics because of the presence of significant instrumental systematics in the 2010 SDO data (private communication by Couvidat, Scherrer, Schou, see also Hathaway et al. 2015).
Appendix B Harmonicanalysis technique
b.1 Harmonic transforms and spectra
Any sufficiently smooth scalar field on the sphere can be expanded in terms of spherical harmonics as
(33) 
where is the orthonormalized spherical harmonic of degree and azimuthal wavenumber , and
(34) 
The onedimensional energy spectrum of is defined as
(35) 
This decomposition can be generalized to vector fields, as explained in the main article; we refer to Eq. (4). The rest of this Appendix describes the procedure used to compute the harmonic transforms of scalar and vector fields derived from SDO/HMI data sets.
b.2 Data interpolation on a GaussLegendreFourier grid
Sphericalharmonics transforms of all fields are computed using GaussLegendre quadrature coupled to a FastFourierTransform algorithm in the azimuthal direction. This requires the knowledge of the field on a twodimensional (2D) GaussLegendreFourier (GLF) grid of size of polar angles corresponding to Gauss nodes , and equallyspaced longitudes (with ; by convention, is the central meridian of the visible disc, is the eastern limb and is the western limb). However, our fields are originally available on a Cartesian 2D grid. Therefore, interpolation from the original Cartesian grid to a GLF grid is a prerequisite. To carry it out, we choose a spectral resolution , compute the Cartesian coordinates of the corresponding GLF grid in the plane of the sky/CCD, and interpolate the field on the original, regular Cartesian grid to these coordinates using the RectBivariateSpline function of the Python module scipy.interpolate. A similar, reverse interpolation procedure is used to reproject fields manipulated in spectral space from the GLF grid to the original “observation” grid. This is, for instance, required to obtain the divergence field on this grid.
b.3 Estimator of spectral coefficients
An apodizing window must be applied to the data in order to deal with the fact that we only see one side of the Sun, and to eliminate nearlimb regions where the CST analysis is unreliable due to projection effects. We use the estimator
(36) 
to obtain the spectral coefficients of a scalar field known on a limited area (see Hathaway (1992) for a similar definition ; the harmonic coefficients of vector fields can be estimated similarly). The normalization coefficients are functionally dependent on the apodizing window, and can be defined in different ways (see next paragraph). The integral in the numerator of Eq. (36), or its generalization to the vector case, are computed from the discrete data sets interpolated on the GLF grid using the python interface of the SHTns library (Schaeffer 2013). The tilde notation for the estimators of the spectral coefficients is dropped in the text to simplify notations.
b.4 Apodizing and normalization issues
Dealing with data on a limited area, or apodizing it, introduces leakage between different harmonics. In helioseismology, this leakage can partly be alleviated by separating the oscillatory components of the signal according to their discrete temporal frequencies (Schou & Brown 1994; Howe & Thompson 1998). However, this approach cannot be used here because the flows that we analyze are characterized by a continuum of time scales. As a result, leakage cannot be eliminated and can notably lead to overestimating the magnitude of the flow.
The results of a detailed technical analysis including careful tests and validations using synthetic data (not presented here but available on request) show that a sensible approach to estimate accurately the spectral coefficients in our problem is to use a “ (heliocentric angle from disc center) window” function,
(37) 
where is a smoothing parameter, is the critical heliocentric angle from disc center above which the apodization becomes significant, and integer parametrizes the sharpness of the falloff of the window function. We adopt an “energyconserving” normalization, which consists in directly calibrating the coefficients introduced in Eq. (36) once and for all based on a comparison between the (known) analytical spectra of a few test synthetic scalar fields defined over the whole sphere on the one hand, and the spectra obtained through a (unrenormalized) harmonic analysis of the apodized versions of the same fields on the other hand. For the family of given by Eq. (37), we found that using a normalization factor independent of and is appropriate. Most importantly, although this strategy does not entirely eliminate the leakage problem, it does alleviate that of overestimating the spectral amplitude.
b.5 Harmonicanalysis parameters
All the results of the paper involving a harmonic analysis were computed for a resolution , corresponding to . This resolution was determined using the Nyquist criterion applied to the CST grid resolution of Mm. We used the window function
(38) 
which apodizes the field sharply at angular distances from disc center larger than , and the same window with a narrower opening . The normalization of the harmonic spectra of the apodized data is based on the energyconserving normalization introduced above, with for all and .