Assessing the Observability of Deep Meridional Flow Cells in the Solar Interior

Meridional circulation regulates the Sun’s interior dynamics and magnetism. While it is well accepted that meridional flows are poleward at the Sun’s surface, helioseismic observations have yet to provide a definitive answer for the depth at which those flows return to the equator, or the number of circulation cells in depth. Here, we explore the observability of multiple circulation cells stacked in radius. Specifically, we examine the seismic signature of several meridional flow profiles by convolving time–distance averaging kernels with mean flows obtained from a suite of 3D hydrodynamic simulations. At mid and high latitudes, we find that weak flow structures in the deep convection zone can be obscured by signals from the much stronger surface flows. This contamination of 1–2 m s−1 is caused by extended side lobes in the averaging kernels, which produce a spurious equatorward signal with flow speeds that are 1 order of magnitude stronger than the original flow speeds in the simulations. At low latitudes, the flows in the deep layers of the simulations are stronger (>2 m s−1) and multiple cells across the convection zone can produce a sufficiently strong signal to survive the convolution process. Now that meridional flows can be measured over two decades of data, the uncertainties arising from convective noise have fallen to a level where they are comparable in magnitude to the systematic biases caused by nonlocal features in the averaging kernels. Hence, these systematic errors are beginning to influence current helioseismic deductions and need broader consideration.


Introduction
The meridional circulation of the Sun is an axisymmetric flow pattern in which material moves from the equator towards the poles at the surface.Because of mass conservation, the flow submerges into the interior, switches direction, and returns towards the equator at some depth.This flow is of great relevance to the global dynamics of the Sun and to its magnetism because it may play an important role in modulating the amplitude and duration of sunspot cycles (e.g., Jiang et al. 2010;Upton & Hathaway 2014a,b) and at transporting angular momentum and magnetic flux (e.g., Wang et al. 1989;Miesch 2005;Featherstone & Miesch 2015).Yet, characterizing the structure of the meridional circulation has been a challenge.
Over the last three decades, many observational attempts have been made to measure the solar meridional circulation (see, Hanasoge 2022, for a recent review).There is a general concensus that at the Sun's surface and in a shallow layer below (i.e., from the photosphere to a depth of ∼ 30 Mm), the meridional flow is predominantly poleward, with a speed of 10-20 m s −1 (e.g., Komm et al. 1993;Švanda et al. 2007;Ulrich 2010;Zhao et al. 2012;Hathaway 2012;Hathaway et al. 2022).Nowadays, improved helioseismic techniques have enabled measurements down to 200 Mm in depth, i.e., near the bottom of the convection zone.However, these extended measurements are inconclusive; different helioseismic analyses of the same data and the same analysis applied to data from different instruments produce different results.In particular, there is evidence for both a single-cell meridional circulation with a deep return flow (e.g., Giles 2000;Rajaguru & Antia 2015;Gizon et al. 2020) and a double-cell meridional circulation with a return flow near the middle of the convection zone (e.g., Hathaway 2012;Zhao et al. 2013;Chen & Zhao 2017;Jackiewicz et al. 2015;Lin & Chou 2018).Interestingly, there is even evidence for both morphologies within the same study.For example, a careful examination of the results of Gizon et al. (2020), suggests single cells during magnetically active time periods and multiple cells during solar minimum (see their Figure 2c).
There are many potential reasons for such inconsistency, some arising from the travel-time measurements and others from the inversion procedure.A prominent issue concerning the measurements is the removal of systematic biases, particularly the so-called center-to-limb effect (Zhao et al. 2012(Zhao et al. , 2016;;Chen & Zhao 2018).Since the physical origin of this bias is unclear, a variety of methods have been suggested for its removal.Complicating this issue is that the systematic biases depend on the mode frequency, radial order, instrumentation, and spectral line choice (e.g.Zhao et al. 2012;Greer et al. 2013;Chen & Zhao 2017, 2018;Rajaguru & Antia 2020).
Many inversion procedures with differing physical assumptions have been developed to measure meridional circulations, and many of the inconsistencies noted previously are likely the result of these differences.For example, some inversion procedures impose mass conservation as a physical constraint on the meridional flow (Schad et al. 2013;Rajaguru & Antia 2015;Mandal et al. 2017;Gizon et al. 2020), and others do not (Zhao et al. 2013;Chen & Zhao 2017).The application of such a constraint does seem to typically result in a single flow cell with depth.However, since direct seismic imaging becomes relatively weak below 0.86  ⊙ , the deep flows in such inversions are largely imposed by an interplay of the constraint and the nature of the inversion's regularization.Since there is a continuum of flow structures that conserve mass, it is unclear whether multiple cells with weak flow could be lurking near the bottom of the convection zone.
In addition to potential assumptions about the flow structure, helioseismic inversions also make assumptions about the nature of acoustic wave propagation.Some procedures employ kernels computed using ray theory, and others use the Born Approximation.The former assumes that the Sun's acoustic waves have infinitely small wavelengths, whereas the latter considers the finite wavelength and scattering of acoustic waves.While, it is indeed true that the averaging kernels produced by the two approximations are different (e.g.Mandal et al. 2018), the two techniques often produce similar results (e.g., see Böning et al. 2017, who used the Born approximation and found that both single and multiple-cell meridional flow profiles can be consistent with their measurements).
Finally, another possibility that has received less attention is that the spatial averaging that is inherent in helioseismic measurements could mask the weaker flows of the meridional circulation in the deep interior.Specifically, the well-known presence of spatial side lobes in the averaging kernels can result in contamination of deep flows with signals from the flows near the surface.In this work, we investigate this possibility by seeking the answer to the following question: If the Sun were to have a meridional circulation similar to those achieved in the numerical simulations, how would the helioseismically observed flow appear?Previous authors have investigated this problem by generating synthetic wavefield data with various meridional-circulation models, and then performed helioseismic measurements to see whether the pre-set flow profiles can be recovered by helioseismic measurements (Hartlep et al. 2013;Stejko et al. 2022).This work uses a complementary approach.
We generate synthetic observations by convolving averaging kernels obtained from the ray approximation (based on those of Zhao et al. 2013) with mean flows obtained from a suite of state-of-the-art numerical simulations of convection in the Sun.We chose this particular set of averaging kernels because they are easily accessible and represent a broad class of currently used inversion techniques.We acknowledge that ray-approximation kernels have limitations (e.g., Birch et al. 2001;Birch & Felder 2004); however, the existence of single versus multiple meridional circulation cells does not depend solely on the choice of ray theory (e.g., Rajaguru & Antia 2015;Böning et al. 2017).Thus, the averaging kernels presented in Zhao et al. (2013) are suitable for exploring the issues we are interested in, such as "side-lobe contamination" and whether or not the deep cell structures that appear in our simulations can be recovered after the convolution.
This paper is organized as follows.In Sections 2 and 3, we describe the meridional flows computed from the numerical simulations and the averaging kernels that we apply to them, respectively.Then, in Section 4 we present and analyze the synthetic observations.Finally, in Section 5 we summarize our results and discuss their implications for our understanding of the deep solar meridional circulation.

Numerical simulations
In this section, we provide a brief overview of our numerical simulations and the resulting flows that we will analyze further in the paper.

Fluid equations
The simulations used in this work correspond to models of the solar convection zone in a 3D spherical shell, including the effects of rotation and density stratification.The simulations were conducted using the open-source Rayleigh convection code (Featherstone et al. 2021), which applies the pseudo-spectral techniques described in Glatzmaier (1984), to solve the anelastic fluid equations (e.g., Ogura & Phillips 1962;Gough 1969)-which are appropriate to study convection and mean flows in the interior of low-mass stars where the flow velocities are substantially subsonic and thermal perturbations are small relative to their mean values.Under the anelastic formulation, the fluid equations can be written in the following form: In the equations above,  is the flow velocity,  = −(  ★ / 2 ) r is the acceleration due to gravity (i.e., we ignore the self-gravity of the convection zone, and assume that the entire mass of the star  ★ lies below the convection zone),  0 = Ω 0 ẑ is the rotation vector (aligned with the axis of the spherical coordinate system, ẑ), D is the viscous stress tensor, () is a function that accounts for radiative heating, and Φ is the rate of viscous heating.Further,  is the mass density,  is the pressure,  is the temperature,  is the specific entropy density, and   is the specific heat capacity at constant pressure.The thermal diffusivity and the kinematic viscosity are denoted by   and , respectively.Note that all thermodynamic quantities are linearized around a temporally steady, spherically symmetric reference state (indicated by an overbar).Fluctuations of these thermodynamic variables about the reference state are denoted by a prime (e.g., the total density is  =  +  ′ ).We approximate the solar convection zone as a perfect gas, so that the reference state obeys  = R , where R is the gas constant.The fluctuations are related by where  is the adiabatic index.
Following Featherstone & Hindman (2016b), we adopt a heating function of the form and the normalization constant  is chosen such that a solar luminosity  ⊙ = 3.846 × 10 33 erg s −1 is deposited in the convection zone For all the models in this work, we adopt impenetrable and stress-free boundary conditions for the velocity field.At the inner boundary, we employ a thermally insulating boundary condition ( ′ / | = inner = 0), and at the outer surface, we enforce that the stellar luminosity  ⊙ deposited into the convection zone by internal heating is carried out via thermal conduction Finally, the physical variables are represented by truncated expansions on the sphere, with spherical harmonics for the angular directions and Chebyshev polynomials for the radial direction.All the simulations are conducted using the same resolution, where the number of radial, latitudinal, and longitudinal points are (  ,   ,   ) = (256, 1536, 3072).After achieving thermal and dynamical equilibration, the models were evolved for at least a thermal diffusion timescale,  2 /  (∼ 10 3 -10 4 rotations in our models).

Nondimensional Control Parameters
Three nondimensional parameters govern the evolution of the flow.These are the flux Rayleigh number Ra, the Ekman number Ek, and the Prandtl number Pr, which in this work are defined respectively as where tildes indicate volume averages over the convective shell.Here,  is the thermal energy flux that derives from the internal heating () (see Featherstone & Hindman 2016b).These numbers can be expressed more concisely in terms of timescales (see e.g., Camisassa & Featherstone 2022), yielding Here,   and   are the viscous and thermal diffusion times across the domain respectively, and  Ω is the rotational timescale.Separately,  ff , which is akin to the convective overturning time, is the characteristic time for a cool parcel to freely fall across the convection zone.
In stellar interiors, these ratios reach extreme values owing to the relative weakness of diffusion when compared to other physical effects.Under solar conditions, order of magnitude estimates give Pr = O (10 −6 ), Ra = O (10 22 ), and Ek = O (10 −15 ) (e.g., Schumacher & Sreenivasan 2020).It is currently impossible, and likely to remain so for decades, to create numerical simulations with solar-like values of the Rayleigh, Ekman, and Prandtl numbers; the spatio-temporal resolution requirements are simply too high.
By contrast, the simulations considered here have Pr = 1, Ra = 3×10 6 , and Ek ∈ [3.6×10 −5 , 2.7×10 −4 ].Fortunately, direct matching of the nondimensional control parameters is not necessary for fruitful study of large-scale flows.Instead, two conditions must be satisfied to ensure that the models and the Sun are in a similar dynamical regime.
First, the Rayleigh number (and so too the degree of turbulence) must be sufficiently large that diffusive processes drop out of the leading-order force and thermal balances.As discussed in Featherstone & Hindman (2016a,b), this condition is satisfied for the Rayleigh and Ekman numbers under consideration here.
The second consideration is that the Coriolis force experienced by convective motions must be sufficiently strong that a solar-like differential rotation, with fast equator and slow poles, is sustained.When this is not the case, a so-called antisolar differential rotation, with rapidly-rotating poles and a slowly-rotating equator will develop.The transition between these two regimes is best characterized by the convective Rossby number, Ro c defined as As long as Ro c is less than unity, a solar-like differential rotation is sustained (e.g., Gastine et al. 2014;Camisassa & Featherstone 2022).Using the solar estimates for our nondimensional control parameters noted above, we find that Ro c = O (10 −1 ).This motivates our choice for the range of Ekman numbers explored in this study, which translates to Ro c ∈ [0.0625, 0.5].
To illustrate qualitatively the differences in the morphology of the flows for different values of Ro c , we first present the two extreme cases in our set of models, namely, Ro c = 0.0625 (strong rotational influence), and Ro c = 0.5 (moderate rotational influence).Figure 1 shows snapshots of the radial velocity near the upper surface (/ outer = 0.968) for both cases.In the low-Ro c case (panel a), the flow structures are noteably thinner and more elongated in latitude than those arising in the high-Ro c case in panel b.In that system, only weak north-south alignment of the flow is apparent, and convective structures possess significantly larger horizontal extent.This variation in convective structure translates into similar variation in the angular momentun transport and, through it, the meridional flow achieved (e.g., Miesch & Hindman 2011;Featherstone & Miesch 2015).In order to analyze the meridional circulation, we examine the mean flows obtained from the simulations.This is done in two steps.First, we average the flows across the azimuthal () direction.Next, we take a temporal average over the periods when the system is in a statistically stationary state (∼ 10 3 -10 4 rotations).Figure 2 shows the latitudinal component of the meridional flow, for the same cases discussed previously (Ro c = 0.0625 in panel a, and Ro c = 1 in panel b).

Flow fields
We observe that the flow consists of multiple cells throughout the shell.However, the key difference between the two runs is that for Ro c = 0.0625, the flows are limited to the surface layers and only at low-latitudes can they penetrate down to approximately half of the convection zone.On the other hand, for Ro c = 0.5, the cells are present throughout the entire depth of the convection zone, with a spherical shape near the surface and a more cylindrical shape towards the interior.
Finally, it is worth mentioning that for the range of Ro c considered in this work, we have not found any solutions where the meridional circulation is composed of a single cell.Although we do find solutions with single cells for Ro c > 1, we exclude them deliberately since they have antisolar differential rotation (e.g., Gastine et al. 2014;Camisassa & Featherstone 2022).
We note that irrespective of the direction of the near-surface meridional flows, which for some of the cases considered here are directed to the poles (resembling the ones in real observations), there exist additional extremely thin circulation cells that hug the outer surface at the equator and at high latitude (a feature that is not present in most observations).These boundary layer cells are confined to the outer 1% of the shell and are the result of the impenetrable boundary condition applied at the upper surface.Such Ekman boundary layers are ubiquitous in rotating systems and should have a thickness proportional to the square root of the Ekman number, Ek.The Sun's Ekman layer is incredibly thin (perhaps 100 m using Ek ∼ 10 −14 ).Our simulations have typical Ekman numbers that are larger than the Sun's, but still small Ek ∼ 10 −4 ; hence, the boundary layers should be 1% of the shell depth, ∼ 2 Mm.When generating the synthetic observations in Section 4, we conduct our analyses with and without the boundary layer.While the Ekman layers in our models are significantly larger than the Sun's, they are still sufficiently shallow that the tiny cells do not change our results.

Averaging Kernels
When using helioseismology to study the meridional circulation in the solar interior, the final product is an estimate of the flow velocity at many positions of interest in latitude and radius (i.e., the target location,   ,   ).We make clear that  represents latitude and not co-latitude.The latitudinal component of this flow can be expressed as where the integral runs over the spatial extent of the convection zone, K is the averaging kernel, and  true is the latitudinal component of the true solar meridional flow.Equation ( 11) can be interpreted as follows: when the flow velocity is measured at a certain location, it represents a spatial average of the actual velocity.Ideally, one would desire that the averaging kernels should be as localized as possible, similar to a  function.However, regardless of the inversion technique used to obtain the flow field, averaging kernels in deeper regions of the convection zone may not be localized, thus complicating the interpretation of the flow field that is obtained.This can result in smoothing over weak features of small size or contamination from distant locations.
We employ averaging kernels that were originally presented by Zhao et al. (2013), who employed time-distance helioseismology (e.g., Duvall et al. 1993) to measure the Sun's meridional circulation throughout the convection zone.The kernels were constructed using the ray approximation, which assumes that the acoustic waves are short wavelength and are sensitive only to the flow field along the raypath that connects two points of observation on the solar surface.Further, only sensitivity to the latitudinal component of the flows is considered, which is justified by the fact that the radial component of the flow has a small influence on the travel-time differences.
Figure 3 shows, for several specific target locations, the averaging kernels that we apply later to our simulated meridional flows.These kernels are well localized near the surface, and become broader as the target depth increases, while still remaining localized in the latitudinal direction1.We emphasize that the averaging kernels possess negative side lobes that appear both above and below the targeted area.This side-lobe behavior is important because when measuring the flow speed deep in the convection zone, this effect introduces contamination from the fast surface flows (i.e., the integral in Equation 11 gets contributions from regions that can be far from the location of interest).This is a major point of this paper, and we investigate this effect in Section 4 by applying the averaging kernels to our hydrodynamical simulations.

Synthetic Observations
We calculate synthetic "observed" flows using Equation (11) with the averaging kernels discussed previously, and assuming that  true is given by the latitudinal component of the mean flows resulting from our simulations, ⟨  ,sim ⟩ (where ⟨•⟩ denotes an average over time and longitude, as explained in Section 2.4).For each of our simulations, the procedure is as follows: 1.For each target location (  ,   ), we make sure that the radial and latitudinal extent of both K and ⟨  ,sim ⟩ cover the same range.Then, since the resolution of the simulation data is larger than that of the kernels, we interpolate the simulated flow fields to the same spatial grid as the kernels using a cubic spline interpolation.
2. Then we compute K (  ,   , , )⟨  ,sim ⟩(, ), and integrate over all  and  as in Equation ( 11).The result is a scalar that corresponds to an estimate of the flow field at (  ,   ).
3.Then, steps 1 and 2 are repeated over the range   / ⊙ ∈ [0.718, 0.981], and As a result, a 2D flow field in the meridional plane is generated.This synthetic observation is the final product of applying the averaging kernels to the simulations.Figure 4 shows a qualitative comparison between the latitudinal component of the meridional flow from the simulations with the synthetic observations.We find that throughout the domain, small-scale features are broadened and, in many cases, disappear due to the spatial averaging imposed by the width of the averaging kernels.We point the reader to the many narrow high-latitude flow features in panels (c) and (d).On the other hand, large-scale flow patterns, particularly in the upper half of the convection zone, are accurately replicated, even when the flow direction reverses near the surface.However, flows in the lower half of the convection zone, even if they survive being smoothed, can be significantly modified by near-surface flows due to side lobes in the averaging kernels.

Discussion
We have confirmed that deep, large-scale flow cells can withstand the convolution process if they are sufficiently fast.This begs the question, "What is the smallest flow speed that can be well-reproduced?".A potential answer can be found when looking into radial profiles of the simulated flows and synthetic observations at different latitudes (see Figure 5, which shows such profiles for Ro c = 0.5, and typifies the results from this analysis for different Ro c as well).
On average, we find that flows whose magnitudes are ≲ 1 m s −1 (largely those at mid and high latitudes in our simulations) are hard to recover after the convolution.They get overwhelmed by a spurious flow of ∼1 m s −1 with the opposite sign.This is clearly seen in the zoomed-in regions of panels b and c.The original flows have speeds of ±0.2 m s −1 and change direction at multiple radial locations.On the contrary, the flows in the synthetic observation have equatorward speeds of roughly −1 m s −1 , and only change their sign in the surface layers.
Only at low latitudes are multiple cells able to survive the convolution process and show up in the synthetic observations (see Figure 4 c-d, and Figure 5 a).There are two reasons for this.First, low-latitude flows are about one order of magnitude faster (∼2 m s −1 ) than flows at the same radii, but higher latitudes.Second, the flows that lie above 0.75 ⊙ , where the kernels have negative side lobes, switch direction with radius and hence tend to cancel out in the convolution process.
Finally, as pointed out by Zhao et al. (2013), the primary difficulty when measuring weak flows are the random errors that pass through the inversion.Estimates of the random errors at low latitude from Zhao et al. (2013) (based on two years of data) are typically ∼1 m s −1 near the surface, and ∼5 m s −1 in the regions below 0.8 ⊙ .We expect that a similar analyses using observations of longer duration would enhance the signal-to-noise ratio throughout the convection zone.For observations spanning a solar cycle, 11 yr, this corresponds to typical errors of ∼1.5 m s −1 .Thus, we find that measurements of the meridional flow averaged over a decade face the unenviable position that the speed and uncertainty of the flow at the base of the convection zone, say 2 ±1.5 m s −1 (e.g.Zhao et al. 2013;Gizon et al. 2020) are both comparable to the systematic bias caused by the negative side lobes in the averaging kernels.We particularly wish to emphasize that the strong poleward flow at the surface (∼20 m s −1 ) results in a spurious equatorward flow at depth (∼1 m s −1 ), potentially overwhelming any additional cells below, complicating the identification of the return flow.While our findings are directly applicable only to the particular inversion procedure in Zhao et al. (2013), they serve as a cautionary note that similar biases may impact other procedures to varying degrees.

Figure 1 .
Figure 1.Radial velocity near the upper boundary (/ outer = 0.968) for the rapidly and slowly rotating cases (panels a and b, respectively).Upflows are indicated in red and downflows in blue.Note that for Ro c = 0.0625, the flow structure varies across latitude, with vortices dominating the poles, and Taylor columns near the equator.On the contrary, for Ro c = 0.5, the convective flows are nearly isotropic.Although the Taylor columns are still present, they are less obvious, and can only be seen near the equator as an alignment of the downflow lanes with the rotation axis.Both cases have the same Rayleigh number Ra = 3 × 10 6 and Prandtl number Pr = 1.

Figure 2 .
Figure 2. Latitudinal component of the meridional flow, averaged over time and longitude as described in the text.Results are shown for the rapidly and slowly rotating cases (panels a and b, respectively).Northward flows are indicated in red and southward in blue.In both cases, there are multiple cells across the shell depth, although for Ro c = 0.0625, the cells are confined to the surface.We zoom in on the flows around the equator to better visualize the tiny cells at the surface which occur in a thin Ekman boundary layer that is typically 2 Mm in thickness.

Figure 3 .
Figure 3. Averaging kernels in Zhao et al. (2013) for three different target positions, located at latitude   = 15 • and at radii   = 0.92 ⊙ , 0.85 ⊙ , and 0.74 ⊙ .For clarity, we only show the low-latitude region of the meridional plane, with  ∈ [ −30 • , +30 • ] and / ⊙ ∈ [0.718, 0.996].Outside this range, the kernels are essentially zero.The complete set of averaging kernels have information for target depths that lie within the range   / ⊙ ∈ [0.637, 0.996].Further the kernels are translationally invariant in the target latitude; meaning that the kernels only depend on the difference between the latitude and target latitude, i.e., K = K (  −   ,   ,  ).

Figure 4 .
Figure 4. Latitudinal component of the meridional circulation in the simulations, ⟨  ,sim ⟩, and synthetic observations,   ,obs .Northward flows are indicated in red, and southward in blue.We make clear that ⟨•⟩ denotes an average over time and longitude.Different panels show results for runs of different Ro c .For Ro c = 0.0625, 0.125, and 0.25, the fastest flows can reach speeds of ∼ ±10 m s −1 , while for Ro c = 0.375 and 0.5, the fastest velocities are ∼ ±20-28 m s −1 .For all the cases presented, there is a good agreement between the simulated and reconstructed flow fields in the surface layers.Differences are more pronounced for deeper flows at mid and high latitudes.

Figure 5 .
Figure 5. Radial profiles of the flow fields for the simulations and synthetic observations.Results are shown for the case Ro c = 0.5.To illustrate the minor effect of the thin boundary layer at the surface (confined to the outer 1% of the shell, we show results for synthetic observations including and removing the boundary layer,   ,obs and  ★  ,obs , respectively) .Panels (a), (b), and (c) show the profiles at  = 15 • ,  = 45 • , and  = 70 • , respectively.We select those positions because they exhibit much weaker flows in the deeper layers of the convection zone (as shown in the zoomed-in regions).