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

The meridional circulation of the Sun is an axisymmetric flow pattern in which material moves from the equator toward the poles at the surface. Because of mass conservation, the flow submerges into the interior, switches direction, and returns toward 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, 2014b) 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 consensus 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⁻¹ (e.g., Komm et al. 1993; Švanda et al. 2007; Ulrich 2010; Hathaway 2012; Zhao et al. 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; Jackiewicz et al. 2015; Chen & Zhao 2017; 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 2(c)).

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, 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 et al. 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 R_⊙, 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 performing helioseismic measurements to see whether the preset 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.

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

2.1. 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:

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot (\overline{\rho }{\boldsymbol{u}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{c}\bar{\rho }\left(\frac{D{\boldsymbol{u}}}{{Dt}}+2{{\boldsymbol{\Omega }}}_{0}\times {\boldsymbol{u}}\right)=-\bar{\rho }{\boldsymbol{\nabla }}\left(\frac{P^{\prime} }{\bar{\rho }}\right)-\frac{\bar{\rho }S^{\prime} }{{c}_{{\rm{p}}}}{\boldsymbol{g}}+{\boldsymbol{\nabla }}\cdot { \mathcal D },\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\begin{array}{l}\overline{\rho }\overline{T}\displaystyle \frac{{DS}^{\prime} }{{Dt}}={\boldsymbol{\nabla }}\cdot \left(\overline{\rho }\overline{T}{\kappa }_{T}{\boldsymbol{\nabla }}S^{\prime} \right)+Q+{\rm{\Phi }}.\end{array}\end{eqnarray} \tag{ 3 }$$

In the equations above, u is the flow velocity, ${\boldsymbol{g}}=-({{GM}}_{\star }/{r}^{2})\hat{{\boldsymbol{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 M_⋆ lies below the convection zone), ${{\boldsymbol{\Omega }}}_{0}={{\rm{\Omega }}}_{0}\hat{{\boldsymbol{z}}}$ is the rotation vector (aligned with the axis of the spherical coordinate system, $\hat{{\boldsymbol{z}}}$), ${ \mathcal D }$ is the viscous stress tensor, Q(r) is a function that accounts for radiative heating, and Φ is the rate of viscous heating. Further, ρ is the mass density, P is the pressure, T is the temperature, S is the specific entropy density, and c_p is the specific heat capacity at constant pressure. The thermal diffusivity and the kinematic viscosity are denoted by κ_T 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 $\rho =\overline{\rho }+\rho ^{\prime}$). We approximate the solar convection zone as a perfect gas, so that the reference state obeys $\overline{P}={ \mathcal R }\overline{\rho }\overline{T}$, where ${ \mathcal R }$ is the gas constant. The fluctuations are related by

$$\begin{eqnarray}&&\frac{\rho ^{\prime} }{\bar{\rho }}=\frac{P^{\prime} }{\bar{P}}-\frac{T^{\prime} }{\bar{T}}=\frac{P^{\prime} }{\gamma \bar{P}}-\frac{S^{\prime} }{{c}_{{\rm{p}}}},\end{eqnarray} \tag{ 4 }$$

where γ is the adiabatic index.

2.2. Model Setup

All simulations in this study use an adiabatically stratified, polytropic reference state designed to resemble the standard solar model of Christensen-Dalsgaard et al. (1996). A detailed description of this reference state is provided in Jones et al. (2011) and Featherstone & Hindman (2016b). We adopt a polytrope of index n = 1.5 (so that γ = 5/3). We set M_⋆ = M_⊙ = 1.989 × 10³³ g, and c_p = 3.5 × 10⁸ erg K⁻¹ g⁻¹. We adopt an inner radius, r_inner, of 0.718R_⊙ and an outer radius, r_outer, of 0.981R_⊙ (where R_⊙ = 6.959 × 10¹⁰ cm is the solar radius). With this choice, the shell depth is L ≈ 183 Mm, the shell aspect ratio is r_inner/r_outer ≈ 0.732 and the mass density varies across the shell by more than two orders of magnitude (ρ_inner/ρ_outer = e⁵ ≈ 148). We set κ_T = ν = 1.175 × 10¹² cm² s⁻¹.

Following Featherstone & Hindman (2016b), we adopt a heating function of the form

$$\begin{eqnarray}&&Q(r)=C\,\overline{P}(r),\end{eqnarray} \tag{ 5 }$$

and the normalization constant C is chosen such that a solar luminosity L_⊙ = 3.846 × 10³³ erg s⁻¹ is deposited in the convection zone,

$$\begin{eqnarray}&&{L}_{\odot }=4\pi {\int }_{{r}_{\mathrm{inner}}}^{{r}_{\mathrm{outer}}}Q(r){r}^{2}{dr}.\end{eqnarray} \tag{ 6 }$$

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 ($\partial S^{\prime} /\partial r{| }_{r={r}_{\mathrm{inner}}}=0$), and at the outer surface, we enforce that the stellar luminosity L_⊙ deposited into the convection zone by internal heating is carried out via thermal conduction

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial S^{\prime} }{\partial r}\right|}_{r={r}_{\mathrm{outer}}}=-{\left.\displaystyle \frac{{L}_{\odot }}{4\pi {r}^{2}\overline{\rho }\overline{T}{\kappa }_{T}}\right|}_{r={r}_{\mathrm{outer}}}.\end{eqnarray} \tag{ 7 }$$

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 (N_r, N_θ , N_ϕ ) = (256, 1536, 3072). After achieving thermal and dynamical equilibration, the models were evolved for at least a thermal diffusion timescale, L²/κ_T (∼10³–10⁴ rotations in our models).

2.3. 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

$$\begin{eqnarray}&&\mathrm{Ra}\equiv \frac{| \tilde{g}| \tilde{F}{L}^{4}}{\tilde{\rho }\tilde{T}{c}_{{\rm{p}}}\nu {\kappa }_{T}^{2}},\quad \mathrm{Ek}\equiv \frac{\nu }{2{{\rm{\Omega }}}_{0}{L}^{2}},\quad \Pr \equiv \frac{\nu }{{\kappa }_{T}},\end{eqnarray} \tag{ 8 }$$

where tildes indicate volume averages over the convective shell. Here, F is the thermal energy flux that derives from the internal heating Q(r) (see Featherstone & Hindman 2016b).

These numbers can be expressed more concisely in terms of timescales (see, e.g., Camisassa & Featherstone 2022), yielding

$$\begin{eqnarray}&&\mathrm{Ra}\equiv \left(\displaystyle \frac{{\tau }_{\nu }}{{\tau }_{\mathrm{ff}}}\right)\left(\displaystyle \frac{{\tau }_{\kappa }}{{\tau }_{\mathrm{ff}}}\right),\quad \mathrm{Ek}\equiv \left(\displaystyle \frac{{\tau }_{{\rm{\Omega }}}}{{\tau }_{\nu }}\right),\quad \Pr \equiv \displaystyle \frac{{\tau }_{\kappa }}{{\tau }_{\nu }}.\end{eqnarray} \tag{ 9 }$$

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 ={ \mathcal O }({10}^{-6})$, $\mathrm{Ra}={ \mathcal O }({10}^{22})$, and $\mathrm{Ek}={ \mathcal 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 spatiotemporal resolution requirements are simply too high.

By contrast, the simulations considered here have Pr = 1, Ra = 3 × 10⁶, and Ek ∈ [3.6 × 10⁻⁵, 2.7 × 10⁻⁴]. 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, 2016b), 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

$$\begin{eqnarray}&&{\mathrm{Ro}}_{{\rm{c}}}={\left(\frac{\mathrm{Ra}\,{\mathrm{Ek}}^{2}}{\Pr }\right)}^{1/2}=\frac{{\tau }_{{\rm{\Omega }}}}{{\tau }_{\mathrm{ff}}}.\end{eqnarray} \tag{ 10 }$$

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 ${\mathrm{Ro}}_{{\rm{c}}}={ \mathcal 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 (r/r_outer = 0.968) for both cases. In the low-Ro_c case (panel (a)), the flow structures are notably thinner and more elongated in latitude than those arising in the high-Ro_c case in panel (b). In that system, only a weak north–south alignment of the flow is apparent, and convective structures possess a significantly larger horizontal extent. This variation in convective structure translates into similar variation in the angular momentum transport and, through it, the meridional flow achieved (e.g., Miesch & Hindman 2011; Featherstone & Miesch 2015).

Zoom In Zoom Out Reset image size

2.4. Flow Fields

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³–10⁴ 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)).

Zoom In Zoom Out Reset image size

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 toward 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⁻¹⁴). Our simulations have typical Ekman numbers that are larger than the Sun's, but still small Ek ∼ 10⁻⁴; 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.

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, θ_T , r_T ). We make clear that θ represents latitude and not co-latitude. The latitudinal component of this flow can be expressed as

$$\begin{eqnarray}&&{u}_{\mathrm{obs}}({\theta }_{T},{r}_{T})={\int }_{\mathrm{CZ}}{ \mathcal K }({\theta }_{T},{r}_{T},\theta ,r){u}_{\mathrm{true}}(\theta ,r){drd}\theta ,\end{eqnarray} \tag{ 11 }$$

where the integral runs over the spatial extent of the convection zone, ${ \mathcal K }$ is the averaging kernel, and u_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 direction. ⁸ 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.

Zoom In Zoom Out Reset image size

We calculate synthetic "observed" flows using Equation (11) with the averaging kernels discussed previously, and assuming that u_true is given by the latitudinal component of the mean flows resulting from our simulations, $\langle {u}_{\theta ,\mathrm{sim}}\rangle$ (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 (θ_T , r_T ), we make sure that the radial and latitudinal extent of both ${ \mathcal K }$ and $\langle {u}_{\theta ,\mathrm{sim}}\rangle$ 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 ${ \mathcal K }({\theta }_{T},{r}_{T},\theta ,r)\langle {u}_{\theta ,\mathrm{sim}}\rangle (\theta ,r)$, and integrate over all r and θ as in Equation (11). The result is a scalar that corresponds to an estimate of the flow field at (θ_T , r_T ).
• 3.  
  Then, steps 1 and 2 are repeated over the range r_T /R_⊙ ∈ [0.718, 0.981], and θ_T ∈ [−75°, +75°]. 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.

Zoom In Zoom Out Reset image size

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).

Zoom In Zoom Out Reset image size

On average, we find that flows whose magnitudes are ≲1 m s⁻¹ (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⁻¹ 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⁻¹ and change direction at multiple radial locations. On the contrary, the flows in the synthetic observation have equatorward speeds of roughly −1 m s⁻¹, 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⁻¹) than flows at the same radii, but higher latitudes. Second, the flows that lie above 0.75R_⊙, 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⁻¹ near the surface, and ∼5 m s⁻¹ in the regions below 0.8R_⊙. We expect that 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⁻¹. 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⁻¹ (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⁻¹) results in a spurious equatorward flow at depth (∼1 m s⁻¹), potentially overwhelming any additional cells below, and 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.

We appreciate the commentary by the reviewers, which led to significant improvements in the manuscript. This work was primarily supported by NASA grants 80NSSC18K1125, 80NSSC19K0267, and 80NSSC20K0193. N.F. and the Rayleigh software were additionally supported by National Science Foundation awards NSF-0949446, NSF-1550901, and NSF-2149126. M.C. acknowledges grant RYC2021-032721-I, funded by the European Union NextGenerationEU/PRTR and AGAUR/Generalitat de Catalunya grant SGR-386/2021. T.H. was supported by NASA contract NAS5-02139. L.K. acknowledges support from NASA through grant No. 80NSSC21K0455. L.I.M. was supported by NSF award AST-2202253. Computations were conducted with support from the NASA High End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center on Pleiades. This work was done in collaboration with the COFFIES DSC, which operates under the cooperative agreement 80NSSC22M0162.