© A. Barekat et al. 2021

1. Introduction

The convection zone (CZ) of the Sun, despite being highly turbulent, shows a well-organized large-scale axisymmetric rotation profile depending on both depth and latitude. The entire CZ rotates faster at the equator than at the poles, and the rotation rate decreases mildly with depth, except near the radial boundaries, where there are regions of strong shear (Thompson et al. 1996; Schou et al. 1998). Additionally, a large-scale circulation in the meridional plane, known as the meridional flow (MC), is also present. The amplitude of the MC is about 15–20 ms⁻¹, which is two orders of magnitude smaller than that of the rotational velocity (Duvall 1979; Hathaway 1996).

The near-surface shear layer (NSSL) occupies about 17% of the CZ, or roughly 35 Mm in depth, from the photosphere. Recently, two further properties of it have been reported. First, the value of the logarithmic radial gradient of the rotation rate is reported to be

$$\begin{array}{c}\hfill \frac{\mathrm{d}ln\mathrm{\Omega }}{\mathrm{d}lnr}\approx -1\end{array}$$(1)

in the upper 13 Mm of the NSSL, independent of latitude up to 60° (Barekat et al. 2014). Second, the gradient evolves over time by an amount between 5–10% of its time-averaged value. This closely follows the magnetic activity cycle (Barekat et al. 2016). On the other hand, the MC maintains its poleward motions throughout the cycle (Hathaway & Upton 2014).

Shear flows play an important role in generating and maintaining the solar magnetic field and its activity cycle (e.g., Krause & Rädler 1980). In particular, radial shear is important in the αΩ dynamo model for explaining the equatorward migration of the magnetic activity (Parker 1955; Yoshimura 1975). In this model, negative radial shear in combination with positive α is required to produce the correct equatorward migration of the activity. This negative shear exists prominently in the NSSL. The effect of the NSSL has been tested numerically in mean-field dynamo models by Käpylä & Korpi (2006), where it was found to aid equatorward migration. More observational and theoretical arguments that the NSSL strongly shapes the solar dynamo process were presented by Brandenburg (2005). The role of NSSL can be easily investigated in mean-field models, where it can be added or removed by hand. In contrast, global 3D convection simulations typically fail to self-consistently generate a realistic NSSL (e.g., Guerrero et al. 2013; Hotta et al. 2015), and thus its role in the resulting dynamo solutions is unclear. Understanding the role that the NSSL plays for the dynamo therefore requires that we understand its formation mechanism and why global simulations do not capture it.

The equations governing the generation of large-scale flows in the solar CZ are the following: First, the azimuthally averaged angular momentum equation describes the time evolution of the differential rotation. This equation is obtained using the Reynolds decomposition, where each physical quantity, A, is decomposed into its mean $\overline{A}$ and fluctuations around the mean, a, and where averages are taken over the azimuthal direction. Then, we obtain the equation

$$\begin{array}{cc}\hfill \frac{\partial }{\partial t}(\overline{\rho }{\varpi }^2\mathrm{\Omega })& =-\mathbf{\nabla }\mathbf{·}\{\varpi [\varpi \overline{\rho {\mathit{U}}^m}\mathrm{\Omega }+\overline{\rho }\overline{u_{\varphi }\mathit{u}}-2\nu \overline{\rho }\overline{\mathit{S}}\mathbf{·}\widehat{\mathit{\varphi }}\hfill \\ \hfill & -({\overline{B}}_{\varphi }\overline{\mathit{B}}/{\mu }_0+\overline{b_{\varphi }\mathit{b}})]\},\hfill \end{array}$$(2)

where $\overline{\rho }$, ${\overline{\mathit{U}}}^m=({\overline{U}}_r,{\overline{U}}_{\theta },0)$, $\mathrm{\Omega }={\overline{U}}_{\varphi }/rsin\theta$, ν, μ₀, and $\overline{\mathit{B}}$ are density, meridional flow, angular velocity, molecular viscosity, the vacuum permeability, and the magnetic field, respectively. Furthermore, ϖ = r sin θ, where θ is the latitude, Q_ϕi and M_ϕi are the Reynolds and Maxwell stresses, and $\overline{\mathit{S}}$ is the mean rate of the strain tensor. The Reynolds and Maxwell stresses are the correlations of fluctuating components $Q_{\varphi i}=\overline{u_{\varphi }{u}_i}$ and $M_{\varphi i}=\overline{b_{\varphi }{b}_i}/{\mu }_0$, respectively, where i denotes r and θ. Density fluctuations are omitted, which corresponds to an anelastic approximation.

Second, the azimuthally averaged equation for the azimuthal component of vorticity describes the time evolution of the MC:

$$\begin{array}{cc}\hfill \frac{\partial {\overline{w}}_{\varphi }}{\partial t}& =\varpi \frac{\partial {\mathrm{\Omega }}^2}{\partial z}+{(\mathbf{\nabla }\overline{s}\times \mathbf{\nabla }\overline{T})}_{\varphi }-{[\mathbf{\nabla }\times \frac{1}{\overline{\rho }}[\mathbf{\nabla }\mathbf{·}(\overline{\rho }\mathit{Q}-2\nu \overline{\rho }\overline{\mathit{S}})]]}_{\varphi }\hfill \\ \hfill & +[\mathbf{\nabla }\times \mathbf{\nabla }\mathbf{·}({\overline{\mathit{B}\mathit{B}}}^T+\mathit{M})],\hfill \end{array}$$(3)

where $\overline{\mathit{w}}=\mathrm{\nabla }\times \overline{\mathit{U}}$ is the vorticity, s is the specific entropy, T is the temperature, and ∂/∂z is the derivative along the rotation axis. The first and second terms describe the centrifugal and baroclinic effects, respectively. From these two equations it becomes clear that meridional flow can drive differential rotation and vice versa, and additionally, any misalignment of the density and temperature gradients can drive meridional circulation through the baroclinic term, while turbulent stresses are important in driving both flows.

Theoretical studies have shown that the main forces generating stellar differential rotation are described by the first two terms in Eqs. (2) and (3) (Rüdiger 1989; Kitchatinov et al. 2013). Additionally, the Coriolis number Ω_⋆, describing the degree of rotational influence on the flow, defined as

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\star }=2\tau \overline{\mathrm{\Omega }},\end{array}$$(4)

where $\overline{\mathrm{\Omega }}$ is the rotation rate of the star and τ is the turnover time of the turbulence, has been found to be a key parameter. Ω_⋆ describes the role of rotation in different parts of the CZ. In particular, it leads to a completely different rotation profile within the NSSL in comparison to the rest of the CZ.

In the solar structure model of Stix (2002), Ω_⋆ changes from the surface to the bottom of the CZ as ${10}^{-3}\lesssim {\mathrm{\Omega }}_{\star }{}^{\text{NSSL}}\lesssim 1\lesssim {\mathrm{\Omega }}_{\star }{}^{\text{CZ}}\lesssim 10$. Helioseismic measurements (Greer et al. 2015) have revealed that in the region of the NSSL, the Rossby number (the inverse of Ω_*) is 0.11 at a depth of 30 Mm and 2.16 in the near-surface layers. This shows that the rotational influence is weak in the surface regions and grows to a significant level in the bottom part of the NSSL. The accuracy of the measurements is generally very good at such small depths.

Nonrotating density-stratified convection is dominated by radial motions, in which case the vertical anisotropy parameter $A_{\text{V}}\propto u_{\text{H}}^2-u_r^2<0$, where $u_{\text{H}}^2=u_{\theta }^2+u_{\varphi }^2$ and $u_r^2$ are the fluctuating horizontal and radial velocities. u_ϕ and u_θ are the longitudinal and latitudinal fluctuating velocities. Rotation tends to suppress convection (e.g., Chandrasekhar 1961), reducing the radial motions, and typically |A_V| decreases when Ω_⋆ increases such that the maximum of |A_V| is achieved for Ω_⋆ = 0 (e.g., Chan 2001; Käpylä et al. 2004). |A_V| is also expected to depend on latitude due to the suppression of $u_r^2$ and the rotational influence on the horizontal components.

Rotation also introduces an anisotropy between the latitudinal and azimuthal directions. This is described by the horizontal anisotropy parameter $A_{\text{H}}\propto u_{\varphi }^2-u_{\theta }^2$. As has been empirically shown using numerical simulations in both local convection (Käpylä et al. 2004) and forced turbulence simulations (Käpylä & Brandenburg 2008), A_H is typically positive and increases with Ω_⋆. Its value also increases toward the equator due the increasingly deviating effect of the Coriolis force on the horizontal flows. Furthermore, A_H → 0 as Ω_⋆ → 0. Thus, Ω_⋆ in the solar CZ also reflects the anisotropy of turbulence, which arises due to the presence of the Coriolis force and density stratification. Consequently, rotation and gravity vectors define the necessary two misaligned preferred directions for nonzero off-diagonal Reynolds stress (Rüdiger 1989).

A theoretical model that reproduces the entire rotation profile of the Sun including the NSSL was presented in Kitchatinov & Rüdiger (2005, hereafter KR05). They derived an extension to the quasi-linear turbulence model of Kichatinov & Rüdiger (1993) that did not take the slow-rotation regime into account that is required to treat the NSSL. In the new model (KR05), the vertical anisotropy parameter A_V is crucial to produce strong inward transport of the angular momentum in the slowly rotating surface layers, resulting in the formation of the NSSL. The horizontal anisotropy parameter is assumed to be small, hence |A_V|≫A_H for Ω_⋆ ≲ 1 applies in the surface layers, although the turbulent transport coefficients depend on latitude. They used a hydrodynamic mean-field (MF) model and studied the value of the anisotropy parameter required to reproduce the NSSL. They concluded that values that agreed reasonably well with the local magnetoconvection models were sufficient to reproduce the NSSL.

The remarkable agreement of the recently observed latitudinal independence of the angular velocity gradient with their model motivated a further development of the theory, in which the effect of the magnetic field in the NSSL was included (Kitchatinov 2016). This led to the prediction that the angular velocity gradient varies during the solar cycle, which qualitatively agrees with the observations. As the variations are caused by the magnetic field, Kitchatinov (2016) suggested that measurements of the rotational properties of the NSSL can be used as an indirect probe to measure the subsurface magnetic field. In their model, however, the Reynolds stresses were computed using a second-order correlation approximation (SOCA), the validity of which in astrophysical regimes with high Reynolds numbers is questionable.

To avoid the necessity of using such simplifications, it is desirable to build numerical simulations of stellar convection, directly solving for the relevant, either hydro- or magnetohydrodynamic, equations in spherical geometry. Such models have been developed and used since the 1970s (e.g., Gilman 1977, 1983; Glatzmaier 1985), but reproducing the NSSL is a serious challenge for these models. These global convection simulations can generate a shear layer close to the equator, mostly confined outside the tangent cylinder, where rotation-aligned, elongated large-scale convection cells form (see, e.g., Robinson & Chan 2001; Käpylä et al. 2011a; Guerrero et al. 2016; Warnecke et al. 2016; Matilsky et al. 2019). Only for a higher-density stratification was a shear layer found that extended to latitudes of 60° (Hotta et al. 2015). In this case, however, the gradient of the radial shear was positive in the range 0° < θ < 45°, which contradicts the helioseismic inferences of the NSSL. Hotta et al. (2015) concluded that the meridional Reynolds stress, originating from the radial gradient of the poleward meridional flow, is the most important driver of the NSSL. In their model, the luminosity was decreased to obtain an accelerated equator, therefore the influence of rotation on convection (Coriolis number) was overestimated, and they also speculated about an unfavorable effect of the boundary conditions on their results. It is therefore unclear whether these results are indeed applicable to the NSSL. Overall, it is cumbersome to use global direct numerical simulations (GDNS) to study the origin of the NSSL because the computational cost is high, many effects are present, and it is difficult to reliably separate them from each other. For such an approach, a simpler modeling strategy is required. This is what we attempt here.

In addition to MF and GDNS models, the NSSL has also been studied from the point of view of different types of equilibria. The most recent of these models, Gunderson & Bhattacharjee (2019), considered the formation of the NSSL in a magnetohydrostatic equilibrium model, where it is driven by a poleward meridional flow near the surface. In addition to the assumption of stationarity, which does not apply because the magnetic field of the Sun is oscillatory, the model considers only nonturbulent states; nevertheless, a large-scale poloidal flow, when inserted in addition to the equilibrium configuration, is seen to reduce the rotational velocity near the surface and therefore leads to an NSSL-like condition there. Miesch & Hindman (2011) accounted for the Reynolds and Maxwell stresses in the governing equations, hence allowing for turbulent effects. The authors considered a case in which an equilibrium condition exists for the angular momentum transport, Eq. (2), in which case the meridional circulation and the relevant stresses must balance. Any imbalance in the term encompassing the stresses was then postulated not only to drive differential rotation, but more importantly, to induce a meridional flow. Similarly, the azimuthal vorticity equation, Eq. (3), in a steady state was postulated not only to drive meridional flow, but more importantly, to contribute to maintaining the differential rotation profile. In the earliest scenario explaining the NSSL, Foukal & Jokipii (1975) proposed that the reason for the existence of the NSSL would be the local angular momentum conservation from rising and falling convective fluid parcels, which would lead to inward angular momentum transport. In the scenario of Miesch & Hindman (2011), however, this angular momentum transport is not a sufficient condition to sustain the NSSL. Another necessary ingredient is the meridional force balance in between the turbulent stresses and centrifugally driven circulation within the NSSL, however. In the bulk of the convection zone, a meridional force balance would be rather provided by the baroclinic effect, and the bottom of the NSSL would be determined by the transition point from baroclinic to Reynolds stress balancing. Some agreement with this scenario was found in the study by Hotta et al. (2015), whose models showed that in the region of the NSSL, the force caused by the turbulent stresses was balanced by the Coriolis force.

Our approach here is entirely different from the approaches reviewed above. We formulate a model with as few ingredients as possible to generate large-scale flows to study the role of rotation-induced Reynolds stress specifically in a rotational regime relevant for the NSSL. This involves replacing convection with anisotropically forced turbulence and omitting density stratification, magnetic fields, and spherical geometry. The simplicity of the model allows an unambiguous identification of the drivers of the mean flows that can be used to assess the generation mechanisms of the solar NSSL. Additionally, this simplicity facilitates measuring the turbulent transport coefficients that can be compared with those in the MF model. This can provide insights regarding the discrepancy between MF and GDNS in obtaining the NSSL.

2. NSSL in terms of mean-field hydrodynamics

In this section we briefly explain the theory of the Λ-effect and its relevance for the formation of the NSSL (Kitchatinov et al. 2013; Kitchatinov 2016). We refer to Rüdiger (1989) for a thorough treatise. In this theory, rotating and anisotropic turbulence contribute to diffusive and nondiffusive transport of angular momentum. The nondiffusive part is known as the Λ-effect (Lebedinski 1941). Therefore the Reynolds stress consists of two parts,

$$\begin{array}{cc}\hfill Q_{\mathit{ij}}& =Q_{\mathit{ij}}^{(\nu )}+Q_{\mathit{ij}}^{(\mathrm{\Lambda })},\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill Q_{\mathit{ij}}& =N_{\mathrm{ijkl}}\frac{\partial {\overline{\mathit{U}}}_k}{\partial x_l}+{\mathrm{\Lambda }}_{\mathit{ijk}}{\mathrm{\Omega }}_k,\hfill \end{array}$$(6)

where N_ijkl and Λ_ijk are fourth- and third-rank tensors describing the turbulent viscosity and Λ-effect, respectively. In spherical geometry, Q_rϕ, Q_θϕ, and Q_rθ are the vertical, horizontal, and meridional stresses, respectively. We note here that the meridional stress appears only in the vorticity equation, and in the model by KR05, it does not play any role in the generation of the NSSL.

Ignoring magnetic fields, the vertical and horizontal stresses are given by

$$\begin{array}{cc}\hfill Q_{r\varphi }& ={\nu }_{\Vert }sin\theta (V\mathrm{\Omega }-r\frac{\partial \mathrm{\Omega }}{\partial r})+{\nu }_{\perp }{\mathrm{\Omega }}^{2}sin{\theta }^{2}cos\theta \frac{\partial \mathrm{\Omega }}{\partial \theta },\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill Q_{\theta \varphi }& ={\nu }_{\Vert }(cos\theta H\mathrm{\Omega }-sin\theta \frac{\partial \mathrm{\Omega }}{\partial \theta })+{\nu }_{\perp }{\mathrm{\Omega }}^{2}sin{\theta }^{2}cos\theta r\frac{\partial \mathrm{\Omega }}{\partial r},\hfill \end{array}$$(8)

where ν_∥ and ν_⊥ are the diagonal and off-diagonal components of the turbulent viscosity tensor N_ijkl, respectively. Component ν_⊥ is caused by the effect of the rotation on the turbulent motions (Rüdiger et al. 2019). V and H are the vertical and horizontal Λ-effect coefficients, which to the lowest order are proportional to A_V and A_H (Rüdiger 1980). These coefficients are typically expanded in latitude in powers of sin²θ as

$$\begin{array}{cc}\hfill V& ={\displaystyle \sum _{i=0}^j{V}^{(i)}{sin}^{2i}\theta ,}\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill H& ={\displaystyle \sum _{i=1}^j{H}^{(i)}{sin}^{2i}\theta .}\hfill \end{array}$$(10)

In the NSSL, Ω_⋆ ≤ 1 and A_H ≈ 0, such that $Q_{\theta \varphi }^{(\mathrm{\Lambda })}$ due to the Λ-effect vanishes. The off-diagonal viscosity ν_⊥ is nonzero, but small, such that its effect is negligible (Rüdiger et al. 2019). It has been shown analytically (Kitchatinov & Rüdiger 2005) and numerically (Käpylä 2019b) that in the slow-rotation regime, only the first term in the expansion of the vertical coefficient V⁽⁰⁾ survives and tends to a constant. Furthermore, when a stress-free boundary condition is applied at the radial boundaries, Q_rϕ = Q_rθ = 0. Using this in Eq. (7) and equating the diffusive and nondiffusive stresses, we obtain

$$\begin{array}{c}\hfill \frac{\partial ln\mathrm{\Omega }}{\partial lnr}=V^{(0)}<0,\end{array}$$(11)

which agrees reasonably well with observational results, where the radial rotational gradient is independent of latitude (Barekat et al. 2014).

3. Model

We investigated the role of Reynolds stresses in the generation of mean flows in a quite simplified model that encompasses all physics required by the MF model to generate the NSSL. The physics are |A_V|≫A_H, Ω_⋆ ≤ 1, and a stress-free boundary condition in the radial direction. The first simplification comes from the geometry, in which we used cubic boxes to map the NSSL at different depths and latitudes instead of using spherical geometry, see Fig. 1. This approach certainly brings specific limitations to this study. They are the absence of the MC and the adoption of stress-free boundary conditions in models that are located at different depths of the NSSL. The consequences of the latter limitation is discussed in detail in Sect. 5.5.

Fig. 1. Schematic representation of the geometry of the current models and their relation to the NSSL. The depth of the layer is magnified for clarity. The simulation boxes are located at nine depths (not all shown) and seven latitudes. Ω⋆ increases gradually from the surface to the bottom of the NSSL.

The second simplification comes from replacing self-consistently generated convection with an artificially forced turbulent medium. We used a forcing that produces anisotropic turbulence, which has two important advantages. First, it allowed us to control the turbulent parameters, for example, the size of the eddies. Second, the desired anisotropy of the NSSL can be provided artificially without gravity. Consequently, the effect of stratification can be simply eliminated from the medium because it is difficult to separate this effect from the effect caused by Reynolds stresses in generating the mean flows.

The third simplification comes from the treatment of thermodynamic properties of the medium. We assumed the medium to be isothermal with a constant sound speed obeying $p=c_s^2\rho$, where p, ρ, and c_s are the pressure, density, and sound speed, respectively. We also note here that full compressibility of the equations is retained, but the effect is weak (Ma = u_rms/c_s = 0.04). This choice leads to a simpler MF analysis in which the density fluctuation can be ignored. This model was used successfully to study the properties of the Reynolds stresses in slow- and fast-rotating regimes by Käpylä & Brandenburg (2008) and Käpylä (2019b).

The hydrodynamic equations that we solved directly are

$$\begin{array}{cc}\hfill \frac{Dln\rho }{\mathit{Dt}}& =-\mathbf{\nabla }\mathbf{·}\mathit{U},\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill \frac{D\mathit{U}}{\mathit{Dt}}& =-c_s^2\mathbf{\nabla }ln\rho -2\mathbf{\Omega }\times \mathit{U}+{\mathit{F}}^{\mathrm{visc}}+{\mathit{F}}^{\mathrm{f}},\hfill \end{array}$$(13)

where D/Dt = ∂/∂t + U⋅∇ is the advective derivative and Ω = Ω₀(cos θ, 0, sin θ)^T is the rotation vector. The viscous force is given by

$$\begin{array}{c}\hfill {\mathit{F}}^{\mathrm{visc}}=\nu ({\mathrm{\nabla }}^2\mathit{U}+\frac{1}{3}\mathbf{\nabla }\mathbf{\nabla }\mathbf{·}\mathit{U}+2\mathit{S}\mathbf{·}\mathbf{\nabla }ln\rho ),\end{array}$$(14)

where ${\mathit{S}}_{\mathit{ij}}=\frac{1}{2}(U_{i,j}+U_{j,i})-\frac{1}{3}{\delta }_{\mathit{ij}}{U}_{k,k}$ is the traceless rate of the strain tensor, δ_ij is the Kronecker delta, and the commas denote differentiation. The forcing function is given by

$$\begin{array}{c}\hfill {\mathit{F}}^{\mathrm{f}}(\mathit{x},t)=\mathrm{Re}(\mathsf{N}·{\mathit{f}}_{k(t)}exp[i\mathit{k}(t)\mathbf{·}\mathit{x}-i\varphi (t)]),\end{array}$$(15)

where x, k, and ϕ are the position, wave vector, and a random phase, respectively. The desired vertical (z) anisotropy can be enforced using a tensorial normalization factor ${\mathrm{N}}_{\mathit{ij}}=({\mathrm{f}}_0{\delta }_{\mathrm{ij}}+{\delta }_{\mathrm{iz}}{cos}^2{\mathrm{\Theta }}_{\mathrm{k}}{\mathrm{f}}_1/{\mathrm{f}}_0){({\mathrm{kc}}_{\mathrm{s}}^3/\delta \mathrm{t})}^{1/2}$ of the forcing, where f₀ and f₁ are the amplitudes of the isotropic and anisotropic parts, respectively. δt and Θ are the time step and the angle between the vertical direction z and k, respectively, and k = |k| determines the dominant size of the eddies. In the forcing, f_k is given by

$$\begin{array}{c}\hfill {\mathit{f}}_k=\frac{\mathit{k}\times \widehat{\mathit{e}}}{\sqrt{{\mathit{k}}^2-{(\mathit{k}\mathbf{·}\widehat{\mathit{e}})}^2}},\end{array}$$(16)

which results in the forcing transversal waves; $\widehat{\mathit{e}}$ is an arbitrary unit vector. The details of the forcing can be found in Brandenburg (2001).

4. Simulation setup

We used the PENCIL CODE¹ (Brandenburg et al. 2021) to run the simulations. We considered a cubic box with size (2π)³ discretized over 144³ grid points. z corresponds to the vertical, x to the latitudinal, and y to the azimuthal direction, respectively, the latter two being referred to as the horizontal directions. Horizontal boundaries are periodic, and stress-free conditions are imposed at vertical boundaries with

$$\begin{array}{c}\hfill U_{x,z}=U_{y,z}=U_z=0\mathrm{on}z=z_{\mathrm{bot}},z_{\mathrm{top}},\end{array}$$(17)

where z_bot and z_top represent the bottom and top of the domain. The box size is represented by the wave number k₁ = 2π/L, and we chose a forcing wave number k_f/k₁ = 10. The units of length, time, and density are $k_1^{-1}$, ${(c_s{k}_1)}^{-1}$, and ρ₀, respectively, where ρ₀ is the initially uniform value of the density. The forcing parameters f₀ = 10⁻⁶ and f₁ = 0.04 were chosen such that the effects of compressibility are weak with a Mach number Ma = u_rms/c_s ≈ 0.04 in all simulations. Moreover, with f₁ ≫ f₀, we fulfilled the NSSL condition in which |A_V|≫A_H; see Table 1. The vigor of turbulence is quantified by the Reynolds number,

$$\begin{array}{c}\hfill \mathrm{Re}=\frac{u_{\mathrm{rms}}}{\nu k_f},\end{array}$$(18)

where $u_{\text{rms}}={(\overline{{\mathit{U}}^2}-{\overline{\mathit{U}}}^2)}^{1/2}$ is the root mean square of the fluctuating velocity field. When a fixed value of the kinematic viscosity is used, $\nu =3.3\times {10}^{-4} {(c_s{k}_1^3)}^{-1}$, the Reynolds number is about 13 for all simulations. We placed the box at seven equidistant latitudes from the equator to the pole by setting the angle θ between the rotation vector and the vertical direction as shown in Fig. 1. The vertical placement was determined by the value of Ω₀, which was varied such that the range of Ω_⋆ from Eq. (4) is relevant for the NSSL. The turnover time was defined as τ = ℓ/u_rms, where ℓ is the size of the eddies. In our simulations the energy-carrying scale of turbulence is the forcing scale ℓ = 2π/k_f. The Coriolis number in the simulations is therefore given by

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\star }=\frac{4\pi {\mathrm{\Omega }}_0}{u_{\mathrm{rms}}{k}_f}.\end{array}$$(19)

The corresponding input parameter is the Taylor number,

$$\begin{array}{c}\hfill \mathrm{Ta}={\left(\frac{2{\mathrm{\Omega }}_0{L}^2}{\nu }\right)}^2.\end{array}$$(20)

The values of Ta, Ω_⋆, and the anisotropy parameters are given in Table 1. An additional run with Ω₀ = 0 was made. This run was used to obtain the reference anisotropy; see Sect. 5.2. It was also used to remove a contribution to the Reynolds stress that appears in the nonrotating case; see Sect. 5.4

Mean quantities were defined as horizontal (xy) averages. The local Cartesian quantities are related to their counterparts in spherical polar coordinates via (r, θ, ϕ)→(z, x, y), $({\overline{\mathit{U}}}_r,{\overline{\mathit{U}}}_{\theta },{\overline{\mathit{U}}}_{\varphi })\to (\overline{{\mathit{U}}_z},\overline{{\mathit{U}}_x},\overline{{\mathit{U}}_{\mathit{y}}})$, Q_θϕ → Q_xy, Q_θr → Q_xz, and Q_rϕ → Q_yz. We normalized the quantities such that ${\stackrel{\sim }{\mathit{U}}}_i={\overline{\mathit{U}}}_i/u_{\text{rms}}$ and ${\stackrel{\sim }{Q}}_{\mathit{ij}}=Q_{\mathit{ij}}/{u_{\text{rms}}}^2$; the tilde denotes this operation. Additionally, the error on the measured physical quantities, which were obtained directly from the simulations, was estimated by dividing the time series into three parts and comparing their time-averaged values with the value obtained from the whole time series. The maximum deviation from the latter was considered to be the error of the measurement.

5. Results

We start our analysis by studying the properties of the velocity field. Then we apply horizontal averages to separate the mean and fluctuating velocities. Using the velocity fluctuations obtained in this way, we measure the diagonal components of the Reynolds stresses to ensure that the flow fulfills the NSSL condition |A_V|≫A_H in Sect. 5.2. In Sect. 5.3, by measuring the mean velocity field at all depths and latitudes, we search for the NSSL-like flow. Next, we focus on the properties of the off-diagonal components of the Reynolds stresses in Sect. 5.4, which are the only transporters of the angular momentum in our model. We also compare their properties with those in previous works of both local and relevant global convection simulations. These comparisons show the differences between our artificially generated turbulent flow with more realistic flows. This can also indicate the relevance of the excluded physics from our model. In Sect. 5.5 we investigate whether these stresses are the main driver of the observed mean flows measured in Sect. 5.3.

5.1. Velocity field

A statistically stationary turbulent state appears after about a few τ independent of Ω_⋆ everywhere except at the equator, where the statistically stationary state is reached between a few to about 300 τ from the lowest to highest Ω_⋆. As an example, we show snapshots of the zonal flow normalized by the sound speed at about 1000 τ for set C46 at the equator and at a latitude of 30° in panels A and B of Fig. 2, respectively. The other components of the velocity field are very similar to the zonal component shown in panel B. The dominant scale of the turbulence is the forcing scale k_f/k₁ = 10. The expected large-scale zonal flow similar to the actual NSSL is generated only at the equator, as shown in panel A. All other sets show a similar behavior.

Fig. 2. Streamlines of the velocity field. The color table shows the amplitude of the azimuthal component of the velocity field normalized by the sound speed. Panels A and B show Uy/cs at the equator and at θ = 30° for set C46, respectively.

5.2. Anisotropy of the flow

We used the diagonal components of the Reynolds stress tensor to measure the anisotropy parameters, which are given by

$$\begin{array}{cc}\hfill A_{\mathrm{V}}& =\frac{Q_{\mathit{xx}}+Q_{\mathit{yy}}-2Q_{\mathit{zz}}}{u_{\mathrm{rms}}^2},\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill A_{\mathrm{H}}& =\frac{Q_{\mathit{yy}}-Q_{\mathit{xx}}}{u_{\mathrm{rms}}^2}.\hfill \end{array}$$(22)

We show the time-averaged diagonal stresses of a nonrotating run in panel A of Fig. 3. All of our runs show similar profiles, regardless of θ and Ω_⋆. The stresses are almost constant in the entire domain, except at the boundaries, where ${\stackrel{\sim }{Q}}_{\mathit{zz}}=0$ and the horizontal components rise to twice higher values. In the interior, the values of ${\stackrel{\sim }{Q}}_{\mathit{zz}}$ are about twice as high as the other two components, reflecting the fact that A_V ≈ −0.5. This is shown in panel B of Fig. 3, where we show volume averages of A_V excluding the boundaries (−2 ≤ zk₁ ≤ 2) to avoid boundary effects. As expected, |A_V| has its maximum value in the absence of rotation at Ω_⋆ = 0.

Fig. 3. Panel A: time-averaged and normalized diagonal components of the Reynolds stress as functions of z. The dotted (solid) line shows ${\stackrel{\sim }{Q}}_{\mathit{zz}}$ (${\stackrel{\sim }{Q}}_{\mathit{xx}}$ and ${\stackrel{\sim }{Q}}_{\mathit{yy}}$) of set C0. The vertical dashed lines mark the part of the domain from which AV and AH were measured. Anisotropy parameters AVpanel B and AHpanel C are shown as functions of Ω⋆ at the latitudes indicated in the legend. The diamonds in panel B show the values of AV at the equator from the runs from which the mean flow was removed.

In the slow-rotation regime, Ω_⋆ ≲ 0.2, the effect of rotation on the turbulence anisotropy is weak at all latitudes, in agreement with the theory. In deeper layers of the NSSL (Ω_⋆ > 0.2), the effect of rotation in generating horizontal motions and suppressing vertical motions becomes more prominent, which is reflected in the two anisotropy parameters. A_H stays positive at all Ω_⋆ and latitudes, which shows that $\overline{u_{\mathit{y}}^2}>\overline{u_x^2}$ independent of the two parameters. We note that the value of $\overline{u_{\mathit{y}}^2}$ increases toward the equator and the bottom of the NSSL. In contrast to $\overline{u_{\mathit{y}}^2}$, $\overline{u_x^2}$ increases toward higher latitudes and at bottom of the NSSL, except at 15°, where it decreases by increasing Ω_⋆. A_H increases significantly by about one order of magnitude from the top to the bottom of the NSSL at low latitudes. In contrast, |A_V| decreases by only about 15% at low latitudes and stays constant at the equator. This shows that the latitudinal dependence of |A_V| is much weaker than that of A_H in Ω_⋆ < 1.

Although the latitudinal dependence of A_V is weak, it shows a nonmonotonous profile. For a specific Ω_⋆, the largest magnitude of A_V is obtained at the equator. |A_V| decreases significantly at a latitude of 15° and then starts to monotonically increase toward higher latitudes, but it never reaches the highest values, which are obtained at the equator. The numbers listed in Table 1 reflect the differences measured between the pole and the equator. They indicate an overall decrease in magnitude of A_V as a function of latitude due to this nonmonotoicity. Käpylä & Brandenburg (2008) applied fully periodic boundary conditions and retrieved a monotonic profile. We performed an additional set of runs at the equator without the mean flows, in which case the monotonic behavior was recovered. The results of these runs are shown in panel B of Fig. 3 and are indicated with a diamond. This shows that the nonmonotonicity is caused by the strong equatorial mean flows generated in our model that are not permitted by the periodic boundary conditions. We note that the GDNS models commonly exhibit thermal Rossby waves (also known as Busse columns or banana cells) close to the equator (e.g., Käpylä et al. 2011b; Guerrero et al. 2013; Hotta et al. 2015; Matilsky et al. 2019), and the mean flows seen in our stress-free boundary condition model might reflect their appearance. Studying this phenomenon is beyond the scope of this study, however and will be addressed at a later time.

Although the anisotropy of the flow changes with increasing rotation rate, we see from comparing the absolute values of A_V and A_H that the model fulfills the NSSL condition |A_V|≫A_H. Therefore the model has all ingredients that were suggested by the MF model as necessary to generate the NSSL.

5.3. Mean flows

The development of mean flows in rotating cases means that it takes significantly longer to reach a statistically steady state than in nonrotating runs. Furthermore, long time averages are needed for a statistical convergence of the turbulent quantities. We ran all the simulations for at least 1100 turnover times. As an example, we show a subset of the time evolution of the three components of the normalized mean velocity field for about 1200 turnover times for set C24 at the equator and at θ = 15° in Fig. 4. At the equator, a large zonal flow $\overline{{\mathit{U}}_{\mathit{y}}}$ with a negative vertical gradient developed gradually over 100τ, as shown in panel A. All other sets show a similar zonal flow profile at the equator, but the amplitude and steepness of the gradient increase with increasing Ω_⋆. Farther away from the equator, the amplitude of the mean zonal flow is significantly reduced and the negative gradient disappears, as shown in panel B of Fig. 4. The dependence of the mean zonal flow on rotation is shown in the panel A of Fig. 5, where we show the time-averaged ${\stackrel{\sim }{\mathit{U}}}_{\mathit{y}}$ at selected Ω_⋆ at a latitude of 15°. When Ω_⋆ increases, the gradient of ${\stackrel{\sim }{\mathit{U}}}_{\mathit{y}}$ changes sign and becomes steeper up to Ω_⋆ = 0.46, then it becomes shallower, and it slowly vanishes in the middle at Ω_⋆ = 1. The latitudinal dependence of ${\stackrel{\sim }{\mathit{U}}}_{\mathit{y}}$ is shown for sets C06 and C46 in panels C and E of Fig. 5, respectively. We find that ${\stackrel{\sim }{\mathit{U}}}_{\mathit{y}}$ decreases as a function of latitude, vanishes at the poles, and the amplitude is lower than 5% of u_rms everywhere except at the boundaries.

Fig. 4. Normalized mean components of the velocity field vs. time in terms of turnover time in representative runs in set C24. The rows from top to bottom show ${\stackrel{\sim }{\mathit{U}}}_{\mathit{y}}$, ${\stackrel{\sim }{\mathit{U}}}_x$, and ${\stackrel{\sim }{\mathit{U}}}_z$. The left and right columns show the mean velocities at the equator and at a latitude of 15°, respectively. To make the comparison of the velocity components feasible, we clipped the values of the color table of panel A at 50% of the maximum value.

The time-averaged meridional component of the mean flow $\overline{{\mathit{U}}_x}$ is consistent with zero at the equator for all runs similar to the run shown in panel C of Fig. 4. In contrast to the zonal flow, its value increases away from the equator; see panel D of Fig. 4. The time-averaged value of this component at 15° is shown in panel B of Fig. 5 for selected values of Ω_⋆. The negative gradient persists up to Ω_⋆ = 0.24. Above this Ω_⋆, the shear slowly vanishes at the center of the box and becomes slightly positive at increasing Ω_⋆. However, the strong shear persists only near the boundaries. We show the latitudinal dependence of $\overline{{\mathit{U}}_x}$ for the two sets C06 and C46 in panels D and F of Fig. 5. The amplitude of $\overline{{\mathit{U}}_x}$ decreases as a function of latitude. The amplitudes of $\overline{{\mathit{U}}_x}$ and $\overline{{\mathit{U}}_{\mathit{y}}}$ are comparable everywhere except at the equator, and the negative gradient of $\overline{{\mathit{U}}_x}$ for Ω_⋆ < 0.1 persists at all latitudes. These results clearly show that in contrast to the MF model, the NSSL-like flow is generated only at the equator. It is also interesting that although our model is significantly simpler than the GDNS models and there is no connection between latitudes, it reproduces results similar to the GDNS.

Fig. 5. Time-averaged normalized mean velocity components vs. vertical direction. panels A and B show ${\stackrel{\sim }{\mathit{U}}}_{\mathit{y}}$ and ${\stackrel{\sim }{\mathit{U}}}_x$ at 15°. The second and third rows show the mean horizontal velocities for sets C06 and C46, respectively, at the latitudes indicated in the legends.

For completeness, we show the vertical component of the normalized mean flow ${\stackrel{\sim }{\mathit{U}}}_z$ in the bottom row of Fig. 4. All runs show a similar pattern of high-frequency oscillations for ${\stackrel{\sim }{\mathit{U}}}_z$ regardless of latitude and Ω_⋆ with amplitudes of about 10⁻⁴u_rms. These oscillations are identified as longitudinal sound waves, as expected for a compressible system in a vertically confined cavity.

5.4. Reynolds stresses

For zero rotation, it is expected that ${\stackrel{\sim }{Q}}_{\mathit{xy}}=0$, see Eq. (8). However, we find that ${\stackrel{\sim }{Q}}_{\mathit{xy}}$ always has a small but nonzero value that also persists in the longest time series of our data. We find that this contribution is present and its magnitude remains unchanged for higher resolutions as well (288³ and 576³ grids). This issue therefore does not appear to be caused by a numerical convergence issue. We have been unable to identify whether the cause is due to some as yet unidentified physical effect, for example, caused by compressibility, the forcing, inhomogeneities in the system, or a combination thereof. Because this contribution is systematically present, we used the nonrotating run (C0) and subtracted ${\stackrel{\sim }{Q}}_{\mathit{xy}}$ from that run from the results of the runs with rotation.

We show representative results of the off-diagonal stresses at five selected Ω_⋆ at 15° (left column) and spatially averaged (−0.5 ≤ zk₁ ≤ 0.5) as function of latitude (right column) in Fig. 6. The Reynolds stresses for all sets are available in Appendix A. The vertical Reynolds stress at all latitudes shows similar profiles as at 15°, see panel A. The stress is nearly constant in the interior of the domain and tends to zero at the boundaries. ${\stackrel{\sim }{Q}}_{\mathit{yz}}$ is always negative, independent of Ω_⋆ and latitude, as shown in panel B. Thus, the vertical angular momentum transport is inward, in agreement with previous studies (e.g., Pulkkinen et al. 1993; Chan 2001; Käpylä et al. 2004; Kitchatinov & Rüdiger 2005; Käpylä & Brandenburg 2008; Käpylä 2019b). Independent of Ω_⋆, the vertical stress vanishes at the pole and has its minimum and maximum amplitude at the equator and 15°, respectively, after which it gradually decreases toward the pole. For a given Ω_⋆, its amplitude is about twice larger at 30° latitude than 60°. The latitudinal dependence of Q_yz is different from that in previous studies by Pulkkinen et al. (1993) and Käpylä et al. (2004) at Ω_⋆ ≈ 1, in which they measured Q_yz from local box convection simulations. Pulkkinen et al. (1993) reported that the latitudinal dependence was almost constant up to 60° and decreased toward higher latitudes. Q_yz has a v-shape profile in latitude with the minimum at 45° in Käpylä et al. (2004). The main ingredient that is missing in our forced turbulence simulation in comparison with theirs is density stratification. Moreover, Käpylä et al. (2004) included the overshooting layer below the CZ. This makes it difficult to determine what causes our results to be different from theirs.

Fig. 6. Left column: time-averaged off-diagonal Reynolds stresses vs. vertical direction at five selected Ω⋆ indicated by the legend at a latitude of 15°. Right column: stresses shown in the left panels, further spatially averaged (−0.5 ≤ zk1 ≤ 0.5), at different latitudes. The rows from top to bottom show ${\stackrel{\sim }{Q}}_{\mathit{yz}}$, ${\stackrel{\sim }{Q}}_{\mathit{xy}}$, and ${\stackrel{\sim }{Q}}_{\mathit{xz}}$.

The middle panels C and D in Fig. 6 show horizontal stress Q_xy. The signature of turbulent fluctuations at the forcing scale is more clearly visible in this component, and the measurement is quite noisy. The values of Q_xy are close to zero up to Ω_⋆ = 0.46, above which they slowly start to become positive (negative) values in the middle (close to the boundaries). We note here that Q_xy obtains notable values at similar Ω_⋆ as A_H. This clearly shows that the generation of horizontal Reynolds stress depends on the presence of the horizontal anisotropy in the flow. However, A_H being maximum at certain latitude does not necessarily mean that Q_xy also has its maximum value at that latitude. This can be seen by comparing the maximum values in panel D of Fig. 6 with the maximum values of A_H shown in panel C of Fig. 3. At a given Ω_⋆, the profile of Q_xy is similar at all latitudes. Its amplitude is maximum at 30° and gradually decreases toward the equator and the pole, as shown in panel D. This result agrees with the observational measurements of Q_xy using sunspot proper motions (Ward 1965; Gilman & Howard 1984; Pulkkinen & Tuominen 1998), but not with the one measured using supergranulation motions, see Fig. 10 in Hanasoge et al. (2016). The horizontal stress has always positive values independent of Ω_⋆ and latitude, in agreement with previous studies in the slow-rotation regime (e.g. Kitchatinov & Rüdiger 2005; Käpylä & Brandenburg 2008; Käpylä 2019b). The latitudinal profile of Q_xy measured by Pulkkinen et al. (1993) is very similar to our results, but has negative values because their box is located in the southern hemisphere; see their Fig. 6.

The meridional stress is shown in the last row of Fig. 6. In contrast to the other stresses, ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ shows a complicated profile, in particular, close to the boundaries. Moreover, it has positive or negative values depending on both Ω_⋆ and θ, as shown in panel E. The latitudinal dependence of the meridional stress is shown in panel F. At Ω_⋆ < 0.1, ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ is positive at low latitudes and ${\stackrel{\sim }{Q}}_{\mathit{xz}}\to 0$ above 45°. By increasing Ω_⋆, ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ moves toward negative values and its absolute value increases. For Ω_⋆ > 0.24, it has its maximum amplitude at about 45° and decreases toward the pole and the equator, similar to ${\stackrel{\sim }{Q}}_{\mathit{xy}}$. The meridional stress in Pulkkinen et al. (1993) also shows a sign change in agreement with ours, but in mid-latitude. However, the sign change occurs at Ω_⋆ ≈ 1 while our results show only negative values at that Ω_⋆. Comparing the absolute amplitude of the stresses in right column of Fig. 6, we see that ${\stackrel{\sim }{Q}}_{\mathit{yz}}$ is always larger than ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ and ${\stackrel{\sim }{Q}}_{\mathit{xy}}$. For example, at Ω_⋆ = 0.64, ${\stackrel{\sim }{Q}}_{\mathit{yz}}$ is about 2 to 10 times larger than ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ and 5 to 20 times larger than ${\stackrel{\sim }{Q}}_{\mathit{xy}}$ depending on the latitudes. When we also compare the absolute amplitude of ${\stackrel{\sim }{Q}}_{\mathit{xy}}$ and ${\stackrel{\sim }{Q}}_{\mathit{xz}}$, we see that $|{\stackrel{\sim }{Q}}_{\mathit{xz}}|>|{\stackrel{\sim }{Q}}_{\mathit{xy}}|$ for all Ω_⋆. These results show that although Q_xy increases as a function of Ω_⋆, its values are still much lower than the vertical stresses, which agrees with the assumption of KR05 regarding the NSSL.

Although our model is quite simple in comparison to the GDNS, it is of interest to compare the Reynolds stresses with simulations such as those in Käpylä et al. (2011b). These authors modeled turbulent convection in a spherical wedge for a variety of rotation rates. Considering the runs of Käpylä et al. (2011b) with Ω_⋆ < 1, we find good agreement for the horizontal stress Q_xy, which is small and positive for small Ω_⋆, and which has appreciable values only for Ω_⋆ > 0.5. However, we find maximu values at 30° instead of at 10…15° in Käpylä et al. (2011b). We also observe a similar trend for Q_xz such that it is positive for small Ω_⋆ in the northern hemisphere with a sign change after certain Ω_⋆. However, this trend depends on latitude in their case; see their Fig. 8. The profile of Q_yz in the convection simulations is quite different from ours, such that it has a strong latitudinal dependence and has both positive and negative values depending on Ω_⋆ and latitude. This is consistent with earlier studies (e.g., Käpylä 2019b), in which a sign change of Q_yz occurs at higher Ω_⋆ than those considered in the present simulations.

5.5. Role of Reynolds stresses in generating the mean flows

As the Reynolds stresses appear in the MF momentum equation, we start by writing the MF equations for $\overline{{\mathit{U}}_x}$ and $\overline{{\mathit{U}}_{\mathit{y}}}$ using Eq. (13). We wrote these equations first by considering that our setup is fully compressible and that the forcing we used here is not solenoidal, which might cause density fluctuations that cannot be ignored in the MF equations. These considerations led to three terms in addition to the Reynolds stresses in the momentum equation (e.g., Käpylä et al. 2020). We compared all of them with the Reynolds stresses. They and their gradients are considerably smaller than the Reynolds stresses. Therefore we can ignore the density fluctuations, and the final set of equations read

$$\begin{array}{cc}\hfill \dot{\overline{{\mathit{U}}_x}}& =-\overline{{\mathit{U}}_z}{\partial }_z\overline{{\mathit{U}}_x}-{\partial }_z{Q}_{\mathit{xz}}-\nu {\partial }_z^2\overline{{\mathit{U}}_x}-2({\mathrm{\Omega }}_y\overline{{\mathit{U}}_z}-{\mathrm{\Omega }}_z\overline{{\mathit{U}}_y}),\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill \dot{\overline{{\mathit{U}}_y}}& =-\overline{{\mathit{U}}_z}{\partial }_z\overline{{\mathit{U}}_y}-{\partial }_z{Q}_{\mathit{yz}}-\nu {\partial }_z^2\overline{{\mathit{U}}_y}-2({\mathrm{\Omega }}_z\overline{{\mathit{U}}_x}-{\mathrm{\Omega }}_x\overline{{\mathit{U}}_z}).\hfill \end{array}$$(24)

Omitting terms proportional to the small quantities ν and $\overline{{\mathit{U}}_z}$, and Ω_y = 0, yields the final form of the equations,

$$\begin{array}{cc}\hfill \dot{\overline{{\mathit{U}}_x}}& =-{\partial }_z{Q}_{\mathit{xz}}+2{\mathrm{\Omega }}_z\overline{{\mathit{U}}_y},\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill \dot{\overline{{\mathit{U}}_y}}& =-{\partial }_z{Q}_{\mathit{yz}}-2{\mathrm{\Omega }}_z\overline{{\mathit{U}}_x}.\hfill \end{array}$$(26)

We verified the validity of the MF equations by considering the steady-state solution, which reads

$$\begin{array}{cc}\hfill \overline{{\mathit{U}}_y}& ={(2{\mathrm{\Omega }}_z)}^{-1}{\partial }_z{Q}_{\mathit{xz}},\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill \overline{{\mathit{U}}_x}& =-{(2{\mathrm{\Omega }}_z)}^{-1}{\partial }_z{Q}_{\mathit{yz}}.\hfill \end{array}$$(28)

We show the horizontal mean velocities in comparison with the RHS of Eqs. (27) and (28) from a latitude of 30° in set C46 in the upper panel of Fig. 7. These results are representative of all nonequatorial cases. Although there are fluctuations in the gradient of the Reynolds stresses, the match is satisfactory.

Fig. 7. Time-averaged mean velocities $\overline{{\mathit{U}}_{\mathit{y}}}$ and $\overline{{\mathit{U}}_x}$ and their corresponding balancing terms in Eqs. (27) and (28) at latitude of 30° (upper panel) and Eq. (30) at the equator (lower panel) in the vertical direction for set C46. In the upper panel, the orange and blue lines show $\overline{{\mathit{U}}_x}$ and $\overline{{\mathit{U}}_{\mathit{y}}}$, respectively. The red and black lines show the RHS of Eqs. (28) and (27), respectively. In the lower panel, the solid and dotted lines show the LHS and RHS of Eq. (30), respectively.

The equator is a particular case, and Eq. (27) cannot be used because Q_xz and Ω_z are both zero there. Therefore we need to use the third component of the MF momentum equation. Applying similar elimination of the terms as for Eqs. (23) and (24), we have

$$\begin{array}{c}\hfill \dot{\overline{{\mathit{U}}_z}}=-c_s^2{\partial }_{z}ln\overline{\rho }-{\partial }_z{Q}_{\mathit{zz}}-2{\mathrm{\Omega }}_x\overline{{\mathit{U}}_y}.\end{array}$$(29)

The pressure gradient appears in this equation due to horizontal averaging. In the steady state, the zonal flow can be written as

$$\begin{array}{c}\hfill \overline{{\mathit{U}}_y}=-{(2{\mathrm{\Omega }}_x)}^{-1}({\partial }_z{Q}_{\mathit{zz}}+c_s^2{\partial }_{z}ln\overline{\rho }).\end{array}$$(30)

We show both sides of Eq. (30) in the lower panel of Fig. 7. The good correspondence indicates that these equations can be used to investigate the role of the stresses in generation of the mean flows.

We emphasize that although in the steady state, for example, at the equator, the two terms on the RHS of Eq. (30) balance, these terms are not the generators of the mean flow. They do, however, determine the final amplitude of the flow. Instead, the mean flows are generated by the gradient of the vertical stress Q_yz at the vertical boundaries, as can be seen from Eq. (26). This flow then slowly penetrates to the middle of domain. This behavior is also clearly shown in the first panel of Fig. 4, where we show the time evolution of $\overline{{\mathit{U}}_{\mathit{y}}}$.

The generation of mean flows is straightforward at the equator because the meridional stress, and hence the meridional flow, vanish there. At other latitudes, the meridional stress and flow have to be included, but it is clear that the Reynolds stresses are the main driver of mean flows in the current setups.

6. Parameterization of Reynolds stresses in terms of mean-field hydrodynamics

Based on the Λ-effect theory explained in Sect. 2, the vertical and horizontal Reynolds stresses given in Eqs. (7) and (8), respectively, can be written in the simulation domain as

$$\begin{array}{cc}\hfill Q_{\mathit{yz}}& =Q_{\mathit{yz}}^{(\nu )}+Q_{\mathit{yz}}^{(\mathrm{\Lambda })}=-{\nu }_{\Vert }\frac{\partial \overline{{\mathit{U}}_y}}{\partial z}+{\nu }_{\Vert }Vsin\theta \mathrm{\Omega },\hfill \end{array}$$(31)

$$\begin{array}{cc}\hfill Q_{\mathit{xy}}& =Q_{\mathit{xy}}^{(\nu )}+Q_{\mathit{xy}}^{(\mathrm{\Lambda })}={\nu }_{\perp }{\mathrm{\Omega }}^{2}sin\theta cos\theta \frac{\partial \overline{{\mathit{U}}_y}}{\partial z}+{\nu }_{\Vert }Hcos\theta \mathrm{\Omega }.\hfill \end{array}$$(32)

Measuring the Λ-effect coefficients V and H from a single experiment is not possible because the turbulent viscosities ν_∥ and ν_⊥ are also unknown. We compensated for this by running another set of otherwise identical simulations, but suppressed the horizontal mean flows artificially at each time step. Therefore the first terms in Eqs. (31) and (32) go to zero. From these simulations, we directly measured Q^(Λ). However, we needed to validate this approach because the velocities can be affected by the nonlinearity of the Navier-Stokes equations. To do this, we performed yet another set of otherwise identical simulations, but used periodic boundary conditions (PBC) in all directions instead of the stress-free boundary condition (SFBC) in the vertical direction. Then we compared the two sets of stresses obtained with these sets of boundary conditions. This comparison of varying boundary conditions is also important with respect to interpreting the Ω_⋆ dependence as a depth dependence. This approach is somewhat artificial because we practically enforced unrealistic BCs within the convection zone.

As an example, we show the horizontal mean velocities for the PBC version of C46 at a latitude of 30° in panel A of Fig. 8. Clearly, no notable mean flow is generated in this run. Therefore the first term in both Eqs. (31) and (32) goes to zero, similar to the cases in which the mean flows are suppressed. In panel B of Fig. 8, we show the results of comparing the Reynolds stresses in PBC and SFBC cases. The difference caused by varying boundary conditions is confined to a very narrow layer near the boundary. These results suggest that our method for the separation of different effects and enforcing artificial SFBC at different depths is valid.

Fig. 8. Panel A: time-averaged normalized mean velocity vs. vertical direction of the PBC run for set C46 at a latitude of 30°. The black and red lines show $\overline{{\mathit{U}}_{\mathit{y}}}$ and $\overline{{\mathit{U}}_x}$, respectively. Panel B: comparison of the time-averaged normalized stresses obtained from periodic and stress-free boundary condition of the same run. The solid and dashed lines show the measured ${\stackrel{\sim }{Q}}_{\mathit{xy}}$ (red), ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ (blue) and ${\stackrel{\sim }{Q}}_{\mathit{yz}}$ (black) from SFBC and PBC runs, respectively.

Considering Eq. (31), the subtraction of the Reynolds stresses obtained from these simulations from the total ones gives

$$\begin{array}{c}\hfill Q_{\mathit{yz}}-Q_{\mathit{yz}}^{(\mathrm{\Lambda })}=-{\nu }_{\Vert }\frac{\partial \overline{{\mathit{U}}_y}}{\partial z}.\end{array}$$(33)

Measuring the vertical gradient of $\overline{{\mathit{U}}_{\mathit{y}}}$, the value of ν_∥ can be determined by performing an error-weighted linear least-squares fit to Eq. (33). Inserting the measured values of ν_∥ back into Q^(Λ) of both Eqs. (31) and (32), we can measure V and H provided that ν_⊥ ≪ ν_∥.

6.1. Properties of the diffusive and nondiffusive parts of Reynolds stresses

Similar to Sect. 5.4, we first measured Q_ij from a nonrotating run and then subtracted its mean value from the corresponding stress in other sets. We show the different contributions to the Reynolds stresses in Fig. 9. In the left column we show stresses from one or two simulation sets, and in the right column we show the dependence of volume averages, excluding the boundaries (−2 ≤ zk₁ ≤ 2), of $Q_{\mathit{ij}}^{(\mathrm{\Lambda })}$ on both latitude and Ω_⋆. In panel A we show the vertical stresses for set C24 at the equator and at a latitude of 30°. With these results, we can explain the minimum of Q_yz at the equator: the diffusive and nondiffusive parts of the stresses are comparable but of opposite signs, leading to a low negative value for the total. With Eq. (31), we see that ν_∥ > 0, which in combination with ${\partial }_z\overline{{\mathit{U}}_{\mathit{y}}}<0$, gives $Q_{\mathit{yz}}^{(\nu )}>0$. Moreover, the final negative value of Q_yz also shows that $Q_{\mathit{yz}}^{(\mathrm{\Lambda })}$ generates the zonal flow. The profile of $Q_{yz}^{(\nu )}$ for all other latitudes is similar to the profile at a latitude of 30°, and it shows that the main contribution from the diffusive part occurs close to the boundaries at |zk₁|≳2 with positive values. Furthermore, the amplitude of $Q_{yz}^{(\nu )}$ decreases toward higher latitudes (not shown). In the middle of the domain, it has negative values, which fits ${\partial }_z\overline{{\mathit{U}}_{\mathit{y}}}>0$ well, as is shown in Fig. 5. The nondiffusive part of the vertical stress is always $Q_{yz}^{(\mathrm{\Lambda })}<0$. Its absolute value increases from the pole toward the equator and increases with Ω_⋆. We also find that $Q_{\mathit{yz}}^{(\mathrm{\Lambda })}$ is linearly dependent on Ω_⋆ in the slow-rotation regime, in agreement with previous numerical results (Käpylä 2019a).

Fig. 9. Panels A, C, and E: time-averaged diffusive and nondiffusive parts of the Reynolds stresses vs. vertical direction. The black and blue lines in panel A show the normalized vertical stresses at the equator and a latitude of 30° for set C24, respectively. In panel C the horizontal stresses are shown at a latitude of 30° for set C64. The blue and black lines in panel E show the meridional stresses for sets C06 and C64 at a latitude of 30°, respectively. The vertical lines denotes the z range used for volume averages. Solid, dotted, and dash-dotted lines show ${\stackrel{\sim }{Q}}_{\mathit{ij}}$, ${\stackrel{\sim }{Q}}_{\mathit{ij}}^{(\mathrm{\Lambda })}$, and ${\stackrel{\sim }{Q}}_{\mathit{ij}}^{(\nu )}$, respectively. Panels B, D, and F: Volume averages over −2 ≤ zk1 ≤ 2, of ${\stackrel{\sim }{Q}}_{\mathit{ij}}^{(\mathrm{\Lambda })}$ vs. Ω⋆ at different latitudes as indicated by the legend.

We show corresponding results for Q_xy in panel C of Fig. 9 for Ω_⋆ = 0.64. $Q_{xy}^{(\mathrm{\Lambda })}$ has positive values in the whole domain, while $Q_{xy}^{(\nu )}$ is almost zero in the middle of the domain, and its contribution to Q_xy is confined to the boundaries at |zk₁|≳2. This also shows that $Q_{xy}^{(\nu )}$ is the main contributor to the negative values of Q_xy close to the boundaries shown in Fig. 6 C. The volume-averaged values of $Q_{xy}^{(\mathrm{\Lambda })}$, excluding the boundaries, are shown in Fig. 9 D as a function of both Ω_⋆ and latitude. Its value is almost zero at the equator and at the pole. It is significantly nonzero above Ω_⋆ > 0.24 and increases with increasing Ω_⋆ independent of latitude. Independent of Ω_⋆, it has maximum value at a latitude of 30°. We note here that the amplitude of $Q_{xy}^{(\mathrm{\Lambda })}$ is also significantly smaller than $\left|Q_{yz}^{(\mathrm{\Lambda })}\right|$. The measured profile of $Q_{xy}^{(\mathrm{\Lambda })}$ is almost identical to the profile obtained by Käpylä (2019b).

Our results for Q_xz are shown in Fig. 9. At low Ω_⋆, $Q_{xz}^{(\mathrm{\Lambda })}$ contributes almost nothing to the total stresses. For Ω_⋆ > 0.15, the contribution of $Q_{xz}^{(\nu )}$ disappears in the middle of the domain, but maintains its positive value close to the boundaries. This is shown in panel E, where we plot Q_xz for low and high Ω_⋆ for the sake of comparison. In panel F, we show volume averages of $Q_{xz}^{(\mathrm{\Lambda })}$ at all Ω_⋆ and latitudes. The value of $Q_{xz}^{(\mathrm{\Lambda })}$ is almost zero at the equator and at the pole. At other latitudes, its absolute value increases with increasing Ω_⋆. It has always negative values, independent of both Ω_⋆ and latitude. These results agree with those of Käpylä (2019b) in the slow-rotation regime.

6.2. Measuring turbulent viscosity

The diagonal turbulent viscosity ν_∥, normalized by ℓu_rms, and its dependence on both Ω_⋆ and latitude is shown in panels A and B of Fig. 10, respectively. Except for the highest latitudes, where measurements are unreliable, the turbulent viscosity decreases monotonically as a function of Ω_⋆ such that for the largest Ω_⋆, corresponding to bottom of the NSSL, its value has decreased by roughly a factor of two. The method used here to measure turbulent viscosity relies on the presence of mean flows. As these diminish toward high latitudes, it is very difficult to obtain reliable estimates of ν_∥ near the pole. We note that the measurements of ${\stackrel{\sim }{\nu }}_{\Vert }$ also suffer from numerical noise at Ω_⋆ < 0.1 at low latitudes. In particular, the latitudinal dependence of ${\stackrel{\sim }{\nu }}_{\Vert }$ for θ ≲ 60°, shown in panel B, is probably not reliable. According to the results at lower latitudes, we conclude that the latitude dependence is weaker than the rotational dependence. As ν_∥ is measured with high confidence, we consider its profile at the equator applicable for other latitudes. We used it to measure V and H at other latitudes. The ratio of turbulent to kinematic viscosity is ν_∥/ν ∼ 10–20, as expected for the fluid Reynolds numbers in the current simulations.

Fig. 10. Normalized turbulent viscosity ${\stackrel{\sim }{\nu }}_{\Vert }$ and Λ-effect coefficients as a function of Ω⋆ and latitude. Panels A, C, and D show ${\stackrel{\sim }{\nu }}_{\Vert }$, V, and H as a function of Ω⋆ from the equator to a latitude of 75°, respectively. Panel B shows ${\stackrel{\sim }{\nu }}_{\Vert }$ as a function of latitude for five selected Ω⋆.

Käpylä et al. (2020) measured turbulent viscosity from nonrotating isotropically forced turbulence with shear. Using the same normalization as in their study with ν_t0 = u_rms/3k_f, we obtain ν_t/ν_t0 ≈ 3.5…3.8 in the slowly rotating simulations in sets C02, C04, and C06. These values are roughly twice higher than those in Käpylä et al. (2020). The cause of the difference is unclear.

We also compared the profile of ν_∥ with an analytical expression for the rotation dependence of the viscosity obtained under SOCA by Kitchatinov et al. (1994, hereafter KPR94). We considered the first term in Eq. (34) of their work, which is relevant to our simulations in which ν_∥ = ν₀ϕ₁(Ω_⋆), where ν₀ = 4/15ℓu_rms is the turbulent viscosity obtained for the isotropic nonrotating case, and where ϕ₁ is a function of Ω_⋆ given in the appendix of KPR94. We scaled the analytical result by a factor of κ = 0.68 to make it comparable with our numerical result. In Fig. 11 we show the result of this comparison. This result shows that except for the κ factor, the rotation dependence agrees fairly well between the theory and numerical simulations.

Fig. 11. Comparison of the obtained turbulent viscosity with the analytical result of KPR94. The solid and dashed lines show the normalized turbulent viscosity and rescaled analytical expression κϕ1, respectively.

Considering the off-diagonal turbulent viscosity ν_⊥, we failed to measure it because the two terms that constitute it, $Q_{xy}^{(\nu )}$ and ${\mathrm{\Omega }}^{2}sin\theta cos\theta \partial \overline{{\mathit{U}}_{\mathit{y}}}/\partial z$, are too small. The measurement error in the former is large.

6.3. Measurements of the vertical Λ-effect coefficient

We measured the vertical Λ-coefficient by substituting the volume averages of $Q_{yz}^{(\mathrm{\Lambda })}$ shown in panel A of Fig. 9 and ν_∥ at the equator using

$$\begin{array}{c}\hfill V=\frac{Q_{\mathit{yz}}^{(\mathrm{\Lambda })}}{{\nu }_{\Vert }sin\theta {\mathrm{\Omega }}_0}.\end{array}$$(34)

Our results are shown in panel C of Fig. 10. |V| is about 0.75 and gradually increases to ≈0.95 for latitudes ≤45°. However, the value of |V| at the lowest Ω_⋆ is lower at all latitudes, but it has large error bars. In contrast to low latitudes, values of |V| at latitudes of 60° and 75° decrease for Ω_⋆ > 0.3. Considering the large errors in the measurements at low Ω_⋆, we might consider V to be roughly constant for Ω_⋆ ≤ 0.15 independent of latitude, but it shows a strong latitudinal and rotational dependence for Ω_⋆ > 0.15. This means that considering the first term V⁽⁰⁾ in Eq. (9) in the NSSL condition alone is not enough, as is assumed in the theoretical model by KR05 explained in Sect. 2. Moreover, the increase in |V| toward higher Ω_⋆ at low latitudes is in contrast with the decrease predicted in KR05 model. The same applies to the results of Käpylä (2019b), who did not consider that ν_t = ν_t(Ω_⋆).

6.4. Measurements of the horizontal Λ-effect coefficient

We measured the horizontal Λ-effect coefficient similarly to the vertical coefficient using

$$\begin{array}{c}\hfill H=\frac{Q_{\mathit{xy}}^{(\mathrm{\Lambda })}}{{\nu }_{\Vert }cos\theta {\mathrm{\Omega }}_0}.\end{array}$$(35)

The results are shown in panel D of Fig. 10. The values of H are always positive, independent of Ω_⋆ and latitude. Its values are one order of magnitude lower than |V| up to Ω_⋆ = 0.6, above which H begins to increase at latitudes < 45°. We also note that its value is zero at the equator and at the pole. H is largest at a latitude of 15° and gradually decreases toward higher latitudes. These results show that close to the surface, H does not play any role in transporting the angular momentum, which validates the assumption applied in the NSSL model by KR05.

7. Conclusions

We applied an alternative approach to the MF and GDNS, that is, we ran direct numerical simulations of forced turbulence in local boxes, to primarily determine whether the assumptions and approximations applied in MF theory to explain the formation of the NSSL are valid. In contrast to the GDNS, we were able to isolate and study the role and contribution of the Reynolds stresses in the rotational regime relevant for the NSSL. Additionally, we were able to measure the turbulent viscosity. Our results show that the three required conditions explained in Sect. 2, which are necessary to generate the NSSL in the RK05 model, are insufficient. In particular, the meridional component of the Reynolds stress cannot be ignored. However, our results are in accordance with Q_xy → 0 in the upper part of the NSSL, whereas Q_xy obtains small but nonzero values close to the bottom of the NSSL, in agreement with theoretical predictions. The role of the vertical Reynolds stress in transporting the angular momentum radially inward agrees with the theory. However, its profile differs from that predicted by the theory. In particular, Kitchatinov et al. (2013) and Kitchatinov (2016) assumed that only the term V⁽⁰⁾ survives in the expansion of V in the NSSL. However, our results indicate that higher-order terms in the expansion of V need to be considered. Moreover, it is also expected from theory that the inward angular momentum flux (V) decreases with increasing Ω_⋆ at all latitudes, but our results show that this expectation is fulfilled only at high latitudes. We also note here that the rotational quenching of the turbulent viscosity, ν_∥, adds another degree of complexity to the problem that was not considered previously in the models of NSSL. From the theoretical MF prediction (Kitchatinov et al. 1994), this behavior agrees qualitatively well with our results, however.

Although these local box simulations have a moderate value of Re ≈ 13 and there is no connection between different latitudes, our results are largely consistent with the stresses and mean flows obtained in the GDNS. On the other hand, the theoretical works used SOCA, which should be valid at Reynolds or Strouhal numbers of up to unity, which is in the vicinity of the parameter regime of the current models. It is expected therefore that we find a relatively good match between the measured turbulent viscosity and the viscosity predicted by SOCA.

Q_xz cannot be disregarded in the NSSL. Its role can be further investigated in more realistic setups using spherical geometry where the artifact of the discontinuity between latitudes can be removed. We also note here that we considered only a single modest Reynolds number and one forcing scale, the effects of which need to be explored with wider parameter studies. The other important physics that needs to be investigated are the effects of stratification, compressibility, and magnetic fields, and this needs to be compared with previous studies that have studied them in turbulent convection, namely Pulkkinen et al. (1993), Chan (2001), and Käpylä et al. (2004).

It is worthwhile to note here that a set of companion laboratory experiments is being proposed to test several aspects of our model. In these experiments, a rotating water-filled apparatus will be used to simulate regions of finite latitudinal extent, including β-plane effects, and forcing will be introduced by pump-driven nozzles at the boundaries (Burin et al. 2019). The relative variation of the system rotation rate and the nozzle exit velocity will allow both the Ω_⋆ > 1 and Ω_⋆ < 1 regimes to be explored. The forcing scale length and isotropy will be changed by opening and closing nozzles and by altering the nozzle shapes and orientations. Time-resolved measurements of the components of the flow velocity will allow the mean flows and stresses to be computed and compared with numerical results and theoretical models. Although the details of the forcing and the fluid boundary conditions will be different in the experiment compared to our simulations, it is expected that meaningful results will be obtained as the rotation rate of the system is varied and the experimental data are analyzed to search for signatures of the Λ-effect.

Acknowledgments

The simulations have been carried out on supercomputers at GWDG and on the Max Planck supercomputer at RZG in Garching. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project “UniSDyn”, grant agreement No. 818665). A. Barekat and E. Gilson acknowledge funding by the Max-Planck/Princeton Center for Plasma Physics. P. J. Käpylä acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG) Heisenberg programme (grant No. KA 4825/2-1).

References

Appendix A: Reynolds stresses

We show the measured ${\stackrel{\sim }{Q}}_{\mathit{yz}}$, ${\stackrel{\sim }{Q}}_{\mathit{xy}}$ and ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ at all latitudes and Ω_⋆ in Fig. A.1, Fig. A.2, and Fig. A.3, respectively.

Fig. A.1. ${\stackrel{\sim }{Q}}_{\mathit{yz}}$ vs. vertical direction at all Ω⋆ indicated by the legend. Each panel shows ${\stackrel{\sim }{Q}}_{\mathit{yz}}$ at a certain latitude.

Fig. A.2. ${\stackrel{\sim }{Q}}_{\mathit{xy}}$ vs. vertical direction at all Ω⋆ indicated by the legend. Each panel shows ${\stackrel{\sim }{Q}}_{\mathit{xy}}$ at a certain latitude.

Fig. A.3. ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ vs. vertical direction at all Ω⋆ indicated by the legend. Each panel shows ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ at a certain latitude.

All Tables

All Figures

	Fig. 1. Schematic representation of the geometry of the current models and their relation to the NSSL. The depth of the layer is magnified for clarity. The simulation boxes are located at nine depths (not all shown) and seven latitudes. Ω⋆ increases gradually from the surface to the bottom of the NSSL.
In the text

	Fig. 2. Streamlines of the velocity field. The color table shows the amplitude of the azimuthal component of the velocity field normalized by the sound speed. Panels A and B show Uy/cs at the equator and at θ = 30° for set C46, respectively.
In the text

Fig. 3. Panel A: time-averaged and normalized diagonal components of the Reynolds stress as functions of z. The dotted (solid) line shows ${\stackrel{\sim }{Q}}_{\mathit{zz}}$ (${\stackrel{\sim }{Q}}_{\mathit{xx}}$ and ${\stackrel{\sim }{Q}}_{\mathit{yy}}$) of set C0. The vertical dashed lines mark the part of the domain from which AV and AH were measured. Anisotropy parameters AVpanel B and AHpanel C are shown as functions of Ω⋆ at the latitudes indicated in the legend. The diamonds in panel B show the values of AV at the equator from the runs from which the mean flow was removed.

In the text

Fig. 4. Normalized mean components of the velocity field vs. time in terms of turnover time in representative runs in set C24. The rows from top to bottom show ${\stackrel{\sim }{\mathit{U}}}_{\mathit{y}}$, ${\stackrel{\sim }{\mathit{U}}}_x$, and ${\stackrel{\sim }{\mathit{U}}}_z$. The left and right columns show the mean velocities at the equator and at a latitude of 15°, respectively. To make the comparison of the velocity components feasible, we clipped the values of the color table of panel A at 50% of the maximum value.

In the text

	Fig. 5. Time-averaged normalized mean velocity components vs. vertical direction. panels A and B show ${\stackrel{\sim }{\mathit{U}}}_{\mathit{y}}$ and ${\stackrel{\sim }{\mathit{U}}}_x$ at 15°. The second and third rows show the mean horizontal velocities for sets C06 and C46, respectively, at the latitudes indicated in the legends.
In the text

Fig. 6. Left column: time-averaged off-diagonal Reynolds stresses vs. vertical direction at five selected Ω⋆ indicated by the legend at a latitude of 15°. Right column: stresses shown in the left panels, further spatially averaged (−0.5 ≤ zk1 ≤ 0.5), at different latitudes. The rows from top to bottom show ${\stackrel{\sim }{Q}}_{\mathit{yz}}$, ${\stackrel{\sim }{Q}}_{\mathit{xy}}$, and ${\stackrel{\sim }{Q}}_{\mathit{xz}}$.

In the text

Fig. 7. Time-averaged mean velocities $\overline{{\mathit{U}}_{\mathit{y}}}$ and $\overline{{\mathit{U}}_x}$ and their corresponding balancing terms in Eqs. (27) and (28) at latitude of 30° (upper panel) and Eq. (30) at the equator (lower panel) in the vertical direction for set C46. In the upper panel, the orange and blue lines show $\overline{{\mathit{U}}_x}$ and $\overline{{\mathit{U}}_{\mathit{y}}}$, respectively. The red and black lines show the RHS of Eqs. (28) and (27), respectively. In the lower panel, the solid and dotted lines show the LHS and RHS of Eq. (30), respectively.

In the text

Fig. 8. Panel A: time-averaged normalized mean velocity vs. vertical direction of the PBC run for set C46 at a latitude of 30°. The black and red lines show $\overline{{\mathit{U}}_{\mathit{y}}}$ and $\overline{{\mathit{U}}_x}$, respectively. Panel B: comparison of the time-averaged normalized stresses obtained from periodic and stress-free boundary condition of the same run. The solid and dashed lines show the measured ${\stackrel{\sim }{Q}}_{\mathit{xy}}$ (red), ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ (blue) and ${\stackrel{\sim }{Q}}_{\mathit{yz}}$ (black) from SFBC and PBC runs, respectively.

In the text

Fig. 9. Panels A, C, and E: time-averaged diffusive and nondiffusive parts of the Reynolds stresses vs. vertical direction. The black and blue lines in panel A show the normalized vertical stresses at the equator and a latitude of 30° for set C24, respectively. In panel C the horizontal stresses are shown at a latitude of 30° for set C64. The blue and black lines in panel E show the meridional stresses for sets C06 and C64 at a latitude of 30°, respectively. The vertical lines denotes the z range used for volume averages. Solid, dotted, and dash-dotted lines show ${\stackrel{\sim }{Q}}_{\mathit{ij}}$, ${\stackrel{\sim }{Q}}_{\mathit{ij}}^{(\mathrm{\Lambda })}$, and ${\stackrel{\sim }{Q}}_{\mathit{ij}}^{(\nu )}$, respectively. Panels B, D, and F: Volume averages over −2 ≤ zk1 ≤ 2, of ${\stackrel{\sim }{Q}}_{\mathit{ij}}^{(\mathrm{\Lambda })}$ vs. Ω⋆ at different latitudes as indicated by the legend.

In the text

	Fig. 10. Normalized turbulent viscosity ${\stackrel{\sim }{\nu }}_{\Vert }$ and Λ-effect coefficients as a function of Ω⋆ and latitude. Panels A, C, and D show ${\stackrel{\sim }{\nu }}_{\Vert }$, V, and H as a function of Ω⋆ from the equator to a latitude of 75°, respectively. Panel B shows ${\stackrel{\sim }{\nu }}_{\Vert }$ as a function of latitude for five selected Ω⋆.
In the text

	Fig. 11. Comparison of the obtained turbulent viscosity with the analytical result of KPR94. The solid and dashed lines show the normalized turbulent viscosity and rescaled analytical expression κϕ1, respectively.
In the text

	Fig. A.1. ${\stackrel{\sim }{Q}}_{\mathit{yz}}$ vs. vertical direction at all Ω⋆ indicated by the legend. Each panel shows ${\stackrel{\sim }{Q}}_{\mathit{yz}}$ at a certain latitude.
In the text

	Fig. A.2. ${\stackrel{\sim }{Q}}_{\mathit{xy}}$ vs. vertical direction at all Ω⋆ indicated by the legend. Each panel shows ${\stackrel{\sim }{Q}}_{\mathit{xy}}$ at a certain latitude.
In the text

	Fig. A.3. ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ vs. vertical direction at all Ω⋆ indicated by the legend. Each panel shows ${\stackrel{\sim }{Q}}_{\mathit{xz}}$ at a certain latitude.
In the text