Generation of zonal flows in convective systems by travelling thermal waves

This work addresses the effect of travelling thermal waves applied at the fluid layer surface, on the formation of global flow structures in 2D and 3D convective systems. For a broad range of Rayleigh numbers ($10^3\leq Ra \leq 10^7$) and thermal wave frequencies ($10^{-4}\leq \Omega \leq 10^{0}$), we investigate flows with and without imposed mean temperature gradients. Our results confirm that the travelling thermal waves can cause zonal flows, i.e. strong mean horizontal flows. We show that the zonal flows in diffusion dominated regimes are driven purely by the Reynolds stresses, always travelling retrograde, while in convection dominated regimes, mean flow advection, caused by tilted convection cells, becomes dominant, which generally leads to prograde mean zonal flows. By means of direct numerical simulations we validate theoretical predictions made for the diffusion dominated regime. Furthermore, we make use of the linear stability analysis and explain the existence of the tilted convection cell mode. Our extensive 3D simulations support the results for 2D flows and thus confirm the relevance of the findings for geopyhsical and astrophysical systems.


Introduction
The problem of the generation of a mean (zonal) flow in a fluid layer due to a moving heat source is an old one. Halley (1687) was probably the first who perceived that the periodic heating of the Earth's surface, due to Earth's rotation, could be the reason for the occurence of zonal winds in the atmosphere. Nearly three centuries later, experiments by Fultz et al. (1959), in which a bunsen flame was rotated around a cylinder filled with water, verified Halley's hypothesis. The moving flame caused zonal flows and the fluid started to move opposite to the direction of the flame. Since then, several experimental and theoretical studies appeared, which illuminated this phenomena.
Thus, Stern (1959) repeated Fultz's experiments using a cylindrical annulus. His observations confirmed the previous results that the fluid acquires a net vertical angular momentum through the rotation of a flame, this time despite the suppression of radial currents in such a domain. Stern then provided a simple two-dimensional (2D) model, showing that the mean motion is maintained through the presence of the Reynolds stresses. Davey (1967) extended Stern's model and provided theoretical explanation that in an enclosed domain, diffusion dominated flows always acquire a net vertical angular momentum in a direction opposite to the rotation of the heat source. His model provided asymptotic scalings for the dependency of the time and space averaged mean horizontal velocity, U x V , with the characteristic frequency of the moving heat source Ω: U x V ∼ Ω 1 for Ω → 0 and U x V ∼ Ω −4 for Ω → ∞. The topic gained further attention when Schubert & Whitehead (1969) suggested that the 4-day retrograde rotation of the Venus atmosphere might be driven by such a periodic thermal forcing. By using a low Prandtl number (P r) fluid, they observed that the induced mean flow rotated rapidly and exceeded the rotation speed of the heat source, which was rotated below a cylindrical annulus filled with mercury (P r 1), by up to 4 times. This validated the linear analysis by Davey, who predicted the speed of the fluid to increase as P r becomes small. However, at this time it became clear that the induced rapid mean flows may exceed the range of validity of Davey's linear theory. Consequently, Whitehead (1972), Young et al. (1972) and Hinch & Schubert (1971) studied the influence of weakly nonlinear contributions. They concluded that the small higher order corrections rather tend to suppress the induced retrograde zonal flows and that the occurring secondary rolls transport momentum in the direction of the moving heat source. It therefore seemed unlikely that the mean flows become much faster than the heat source phase speed, even for small P r, as soon as convective processes come into play.
The preceding analysis certainly lacked the complexity of convective flows, and therefore Malkus (1970), Davey (1967) and other authors anticipated that convective and shear instabilities could alter the entire character of the solution. In particular, Thompson (1970) showed that the interaction of a mean shear with convection can lead to a tilt of the convection rolls and thus to the transport of the momentum along the shear gradient and thereby amplifies the mean shear flow. In this scenario, the convective flow is unstable to the mean zonal flow even in the absence of a modulated travelling temperature variation, which suggests that the mean zonal flows might be the rule and not the exception to periodic flows that are thermally or mechanically driven. However, the direction of this mean zonal flow would be solely determined by a spontaneous break of symmetry; it could either move counter (retrograde) to the imposed travelling wave (TW) or in the same directions as the TW (prograde).
The existence of mean flow instabilities in internally heated convection and in rotating Rayleigh-Bénard convection (RBC) (Ahlers et al. 2009) was studied theoretically by Busse (1972Busse ( , 1983 and Howard & Krishnamurti (1986), but has not been observed in laboratory experiments. In classical RBC, a zonal flow, if imposed as an initial flow, can survive (Goluskin et al. 2014) but only if the ratio of the horizontal to vertical extensions of the domain is smaller than a certain value, see Wang et al. (2020b) and Wang et al. (2020a). Also, several studies examined the effects of time-dependent sinusoidal perturbations in RBC. Venezian (1969) showed that the onset of convection can be advanced or delayed by modulation, while Yang et al. (2020) and Niemela & Sreenivasan (2008) demonstrated a strong increase of the global transport properties in some cases.
Its general nature makes the travelling thermal wave problem appealing to study, however, to our knowledge, there are only a few studies recently published, that are related to the original "moving flame" problem. Therefore in the present study we revisit the existing theoretical models, specifically Davey's model and validate it by means of state of the art direct numerical simulations (DNS). Furthermore, we study a setup with non-vanishing vertical mean temperature gradient (as in RBC), to study the influence of the travelling thermal wave on convection dominated flows and discuss the absolute strength and the direction of the induced zonal flows. Despite the substantial advances over the years, it remains unanswered, whether the thermal travelling wave problem is merely of academical interest or, indeed, of practical relevance in the generation of geoand astrophysical zonal flows (Maximenko et al. 2005;Nadiga 2006;Yano et al. 2003).
For this purpose, in chapter 2, we complement our analysis with thorough 3D DNS. For the sake of generality, we choose a classical RBC setup. Ultimately, we analyse the absolute angular momentum in 3D flows (respectively, horizontal velocity in 2D flows) and provide insight into the mean flow structures.
Here t denotes time and e z the unit vector in the vertical direction. The equations have been non-dimensionalised using the free-fall velocity u f f ≡ (αg∆Ĥ) 1/2 , the free-fall time t f f ≡Ĥ/u f f , ∆ the amplitude of the thermal TW andĤ the cell height. The dimensionless parameters Ra, P r and the aspect ratio Γ are defined by: whereL is the length of the domain, ν is the kinematic viscosity, α the isobaric thermal expansion coefficient, κ the thermal diffusivity and g the acceleration due to gravity. This set of equations is solved using the finite-volume code goldfish (Kooij et al. 2018;Shishkina et al. 2015), which employs a fourth-order discretization scheme in space and a third order Runge-Kutta, or, alternatively, a Euler-Leapfrog scheme in time. The code runs on rectangular and cylindrical domains and has been advanced for a 2D-pencil decomposition for a highly parallel usage. In this study the following notations are used: Temporal averages are indicated by an overline or by a capital letter, thus the Reynolds decomposition of the velocity reads u = U +u , decomposing u into its mean part U and fluctuation part u . Unless specifically stated, time averages are carried out over a long period of time, however, in section 3.1.1, the averaging period was deliberately restricted to only a few wave periods to achieve a time scale separation. Further, the spatial averages are denoted by angular brackets · , followed by the respective direction of the average, e.g. · x denotes an average in x; · V denotes a volume average. And ultimately, the velocity vector definitions u ≡ (u x , u y , u z ) ≡ (u, v, w) are used interchangeably.

Theoretical model
Already the earliest models proposed by Stern (1959) and Davey (1967) gave a considerable good understanding of the moving heat source problem. Although there are more complex models (Stern 1971) based on adding higher order non-linear contributions (Whitehead 1972;Young et al. 1972;Hinch & Schubert 1971;Busse 1972), this section focuses on revisiting the main arguments of Davey's original work, which is expected to give reasonably good results in the limit of small Ra. Besides, a more complete derivation and concrete analytical solutions are provided in Appendix A.
Given the linearized Navier-Stokes equations in two dimensions and averaging the horizontal momentum equation in the periodic x-direction and over time t, one can derive the following balance: Evidently, a mean zonal flow U x is maintained by the momentum transport due to the Reynolds stress component u w and by mean advection through W ∂ z U . The theory further advances by assuming that no vertical mean flow exists (W = 0), which reduces equation (2.1) to the balance between viscous mean diffusion and Reynolds stress diffusion. Furthermore, by neglecting convection and variations in x, the linearized equations can be written as a set of ordinary differential equations, that can be solved sequentially to find u and w and ultimately the Reynolds stress term u w . This procedure is shown in Appendix A. Given the Reynolds stress field, equation (2.1) has to be integrated twice to obtain the mean zonal flow U (z). Integrating that profile again finally gives the total mean zonal flow U V , which is an important measure of the amount of horizontal momentum or, respectively, angular momentum in cylindrical systems, that is generated due to the moving heat source. The last step can be solved numerically, however, following Davey (1967), the limiting relations can be calculated explicitly: where the horizontally travelling wave, θ(x, t) = 0.5 cos(kx − 2πΩt), is applied to the bottom and top plate. Here k denotes the wave number of the TW. We would like to add that this theoretical model is, as determined by its assumptions, expected to be limited to diffusion dominated, small-Ra flows. However, when momentum and thermal advection take over, its validity remains questionable. We will show later that after the onset of convection, where eventually mean advection takes over, the neglect of the W ∂ z U -contribution is no longer justified.

2D-convective system
As described by Stern (1959), the generation of a laminar zonal flow by a TW can be successfully explained in a 2D system, which makes it a good starting point. The Here, Ω indicates the temporal frequency of the travelling TW in free-fall time units. For example, Ω = 10 −1 describes a wave with a period of 10 free-fall time units τ f f , and ∆θ is introduced as a control parameter for the strength of the mean temperature gradient. In the following, two different setups are considered. In setup A (figure 1 a) -the one originally examined by Davey (1967) -no mean temperature gradient exists (∆θ = 0) and the top and bottom plate temperatures are equal, whereas in setup B (figure 1 b) a mean, unstable temperature gradient is applied (∆θ = 1). For simplicity, the mean temperature gradient is set equal to the amplitude of the thermal wave. In this setup, effects of convection are expected to become dominant. Averaged over time, this setup resembles RBC, therefore, it can be regarded as a spatially and temporally modulated variant of RBC. Further, no-slip conditions are applied at the top and bottom plates, the x-direction is periodic and the domain has the length L = 2π and the height H = 1.
The overall focus in this study lies on variations of the zonal flow with Ra and Ω. Thus, the parameter space spans 10 3 Ra 10 7 and 10 −4 Ω 10 0 , while the aspect ratio and Prandtl number are kept constant (Γ = 2π, P r = 1). Exemplary temperature fields at a fixed Ω = 0.1 are shown in figure 2.

Results
The theoretical model, as presented in Appendix A, aims to explain the generation of the total mean momentum U x V for a given Ra and wave frequency Ω. Moreover, it predicts that the generated mean momentum will be directed opposite, i.e. retrograde, to the travelling thermal wave. In this section we study the validity of the model and reveal its limitations.  Figure 3. Mean velocity of the zonal flow vs. the wave frequency Ω for Ra = 10 3 ( ), 10 4 ( ), 10 5 ( ), 10 6 ( ) and 10 7 ( ). Circles (stars) denote a retrograde (prograde) mean zonal flow, the solid lines of the corresponding colour show the results of the theoretical model by Davey (1967). (a) Setup A and (b) Setup B. Figure 3 shows the numerical data from the DNS together with the respective results of the theoretical model, for different Ra. Worth noting first is, that the maximum of the theoretical model is located at a fixed frequency, if the frequency is expressed in terms of the diffusive timescale rather than the free-fall time scale Ω κ,max = Ω/ √ RaP r ≈ 0.66. This indicates, that the model predictions could be collapsed onto a single curve. Nonetheless, this was avoided here for the sake of clarity.
We begin our discussion with the results of setup A, shown in figure 3 (a). The theoretical model by Davey (1967), indicated by the solid lines, gives accurate results for Ra = 10 3 and a good agreement for Ra = 10 4 , although, evidently the model systematically overestimates the mean momentum generation for higher Ra. In fact, this is consistent with Whitehead (1972), Young et al. (1972) or Hinch & Schubert (1971) who observed that corrections of higher order non-linear contributions tend to suppress the induced retrograde zonal flows. Also it suggests that an induced mean flow does not strengthen itself, i.e. there is no positive feedback mechanism between the mean flow and Reynolds stresses. While all low Ra flows and high Ra flows in the limit of large Ω are well predicted by the model, the large Ra flows are mostly over predicted (except Ra = 10 7 /Ω = 0.1, the only flow of that setup that becomes truly turbulent, despite similar initial conditions). Presumably even more important is that some of the flows for Ra 10 5 exhibit a positive/prograde mean flow, indicated by a star symbol, which is especially prevalent at small Ω.
Turning the focus to Setup B, shown in figure 3 (b), the differences become even more obvious, since adding a mean temperature gradient enhances the effects of convection further. For Ra = 10 3 the picture is clear, as it is below the onset of convection Ra c ≈ 1708 for classical RBC even for the unbounded domains. The Reynolds stresses remain dominant, which preserves the development of a mean flow opposite to the TW direction. However, for Ra 10 3 all but a few of the simulations end up with a prograde mean flow final state. In order to understand the role of the mean flow, we analyse the two terms on the right side of equation (2.1), which are presented in figure 4. The model neglects mean advection, it only captures contributions of u w . As seen in figure 4 (a), this is justified for a flow without strong convection effects and the model predictions agree well with the Reynolds stresses obtained in the simulations. This is different from the situation in figure 4 (b), where obviously mean flow advection W ∂ z U starts to take over. The shape of the mean flow advection curve is antiphase to the Reynolds stress curve The underlying reason for that will be examined in more detail in the next section. But briefly, the main argument is that there exist two competing mechanisms, one induced by the TW and the other induced by convection rolls, which act on different time scales. At small Ra, as convection rolls move considerably slower, an average over a few TW time periods can reliably separate both structures, so that the Reynolds stresses reflect mainly the TW contributions, while the mean field represents the convection rolls. Therefore, the dominant mean flow advection in figure 4 (b) reflects the dominance of advection by convection rolls as Ra increases.
A few more interesting observations can be deduced from figure 3 (b). First, compared to the theory, the simulations show significantly larger values at high Ra. Apparently, the mean zonal flow can be substantially stronger than expected and its velocity can exceed the TW phase velocity. Second, while the theory predicts the location of the maximum zonal flow at a constant diffusive time scale, the DNS indicates a coupling with the freefall time rather than with the diffusive time and the maximum is found in the region 0.01 Ω 0.1. This is in so far important, since natural flows often fall within this parameter range. We show this in the context of the Earth's atmosphere in section 4. Finally the instantaneous fields (figure 3 b) most often show three plumes (figure 2 b), while a classical RBC simulation with the same initial conditions would develop four plumes. Presumably, either the sinusoidal temperature distribution at the plates, or a preexisting shear flow (before Rayleigh-Taylor instabilities develop) reduces the number of plumes. On this basis, we tested the linear stability of the Rayleigh-Taylor instabilities with an imposed shear flow, and found indeed that the wavelength of the most unstable mode decreases.

Origin of prograde flows in convection dominated flows
In order to understand how prograde flows can emerge, we looked at the route from small to large Ra for a specific configuration. Setup B and Ω = 0.1 is well suited for this purpose, since the transition from a retrograde flow to a prograde flow appears early, already below Ra = 10 4 (figure 3 b). Thus, a simulation was initiated at Ra = 1000 and then Ra was progressively increased by 1000 each time after a steady state has settled. The time evolution of the total mean zonal flow is given in figure 5 (a). At the lowest Ra, the mean flow is retrograde. Increasing Ra to 2000 enhances its strength further, as anticipated. But already at Ra = 3000 the zonal flow breaks down and its vertical profile, as seen in figure 5 (b), flattens. Ultimately, at Ra 4000 this profile flips over and the total zonal flow turns into a prograde state.
As we have shown in the preceding analysis (figure 4 b), in the presence of convection cells, the mean zonal flow can be fed by the base flow itself, in particular it is fed by the vertical advection of horizontal momentum W ∂ z U . Now, let us consider perfectly symmetric convection cells; although locally, at a position in x, momentum may be transported up-or downward, the symmetry, however, would balance this transport at another location and the net transport would become zero. Therefore there must be a symmetry breaking in the convection cells, which correlates W with ∂ z U . A possible mechanism, even discussed in the context of the moving heat source problem, was described by Thompson (1970) and theoretically analysed by Busse (1972), who showed that, in a periodic domain, convection rolls can become unstable to a mean shear flow. This mean shear tilts the convection cells such that their asymmetric circulation maintains a shear flow. In the following this mean flow instability will be called tilted cell instability. Busse (1972) showed the existence of this instability on a analytic base flow field. Differently, in the following we conduct a stability analysis on a base flow extracted from the DNS.
The first rise of the curve in figure 5 a at Ra = 3000 coincides with the observed onset of convection, which is slightly delayed compared to classical, unmodulated RBC (Ra c ≈ 1708). The convection cells at that point appear to be standing still, almost unaffected from the TW and clearly orders of magnitudes slower than the TW. Therefore a short time average, over one wave period, was applied to separate both timescales, which results in the base flow as shown in figure 5 c. Based on this base flow, a linear, temporal stability analysis of the full 2D linearized Navier-Stokes equations was conducted. Details herefore are given in the Appendix C. While no unstable mode was detected for Ra = 3000, for Ra = 4000 the mean flow becomes unstable, to the mode presented in figure 5 d. The growth rate of it is σ ≈ 0 + 0.2i, suggesting no oscillatory behaviour (real part is zero) but exponential temporal growth (imaginary part larger than zero). This mode shares characteristics with the tilted cell instability described by Thompson (1970), in the sense that the mode induces a mean shear flow (see profile in figure 5 d). However, rather than the "pure" shear flows as presented by Thompson (1970) and Busse (1972Busse ( , 1983) with a vanishing total net momentum when integrated vertically, the fluctuation profile found in our study (figure 5 d on the right) shows a more directed flow, negative in the vicinity of the plates and stronger positive in the center. And especially interesting, its momentum profile has a similar shape as the final state solution of typical prograde flows, e.g. the profile on the right in figure 5 b. A few more notes are necessary. The difference between the shape of the mode found in this work, compared to the ones from Thompson and Busse might be explained by different BCs, as both authors applied free-slip conditions at the plates, in contrast to our no-slip conditions. In addition, in their seminal works and in the work of Krishnamurti & Howard (1981), it was already remarked that the mean flow transition is caused by a spontaneous symmetry break and therefore the direction of the shear flow is somewhat arbitrary as it depends on the initial conditions. Indeed, a change in the grid size of the stability analysis led to a most unstable mode with a reversed shear flow profile compared to the mode shown in figure 5 d. And finally, even though in figure 5 a tilted rolls are shown to start later as convection rolls, it actually is likely that the convection cells tilt as soon as convection sets in, it is just not clearly visible from the flow fields at that point.
In a nutshell, the mean flow is unstable -even in the absence of a boundary temperature modulation -to a mode with tilted convection cells and non-zero total mean horizontal velocity. Both modes, prograde and retrograde, are found in the global stability analysis, thus it remains unanswered why the DNS at high Ra almost exclusively end up moving in the same direction as the TW. The disturbance velocity profiles resemble those of the final mean flow velocity profiles, therefore, the presented mean flow instability is a plausible mechanism for the generation of moderate strong zonal flows after onset of convection, then dominating over the Reynolds stress mechanism, that is inherent to diffusion dominated flows.

Space-time structures
The flows found in this study revealed surprisingly rich formations. Therefore this part will be completed with examples of some space-time structures that have been observed in the 2D system and, already ahead of the next part, in the 3D cylindrical system. In addition, movies are provided as supplementary material.
In general, in 2D, as can be seen from figure 2 the temperature field is either symmetric around the horizontal mid-plane (Setup A), or not; in this case there exist plumes (Setup B). In the latter case, there are usually three up-and three down-welling plumes identifiable. In the 3D case, the flow consists of rising and falling plumes, which together form a large scale circulation (LSC). If the TW propagates slowly (small Ω), the plumes (2D) or respectively the LSC plane (3D) drift with the same speed as the TW and both structures appear to be connected. However, as Ω increases and, hypothetically, the TW timescale τ Ω becomes small compared to thermal diffusion τ κ (τ κ /τ Ω = √ P rRaΩ > crit), the plumes (2D) or LSC (3D) "break-off" from the TW, forming two separate structures, acting on different timescales. Figure 6 shows the space-time structures of the temperature field, evaluated at midheight, and in the 3D case at mid-height and near the sidewall. The structures at midheight either (i) travel with the same speed (but a phase difference) as the thermal wave (a,d), or travel with phase speeds different to the thermal wave and in this case either (ii) retrograde (b,e) or (iii) prograde (c). Regime (i) is expected for small Ra and/or small Ω parameters, (ii) is found for large Ra and large Ω, if no mean temperature is present, and (iii) exist in strongly convection dominated flows for large Ra and large Ω, especially if a mean temperature gradient is present. Furthermore it is striking that temperatures between the left and the right region in the vicinity of the plumes center (hottest or coldest regions in figure 6) do not necessarily fill with the same temperature (c). This gives further evidence for a mean flow instability, as it features similarities of the temperature field of the unstable mode given in figure 5 (d), due to which a plume loses its horizontal symmetry. Considering the speed of the drifting plumes (b,c), we observe initially exponential growth, as anticipated from an instability, followed by a, possibly, non-linear saturation.

3D-convective systems
The preceding part, as most of the existing literature, is confined to 2D flows. Now we will discuss the moving heat source problem in the context of more complicated 3D convective flows. In general, travelling wave solutions are common amongst 3D convective systems. Bensimon et al. (1990); ;  observed convection rolls propagating azimuthally in a large aspect ratio annulus near the onset of convection. Their drift velocity is of the order of magnitude 10 −4 to 10 −3 , however, drift velocity is not necessarily equal to the mean azimuthal flow. Another kind of travelling wave solution in RBC systems are the spiral patterns found in large aspect ratio cells (Bodenschatz et al. 1991(Bodenschatz et al. , 2000. These spirals are rotating in either direction, although corotating spirals are more numerous (Cross & Tu 1995), and are known to be coupled with an azimuthal mean flow (Decker et al. 1994). Furthermore, in rotating systems travelling wave structures are quite common (Knobloch & Silber 1990). These structures are strongly geometry-dependent (Wang et al. 2012) and known to induce mean zonal flows that propagate pro-and retrograde (Zhang et al. 2020).
Despite the vast literature on these phenomena, quantitative data on mean flows that are induced by external travelling thermal waves in 3D flows seems to be rare. Therefore our main goal in this part is to gain insight on the strength and structure of such mean flows, and discuss whether their order of magnitude is relevant in natural flows. For this purpose we took the paradigm convective system cylindrical RBC and studied it by means of direct numerical simulations.

Numerical Setup: Cylindrical Rayleigh-Bénard Convection (Pr=1)
The setup is essentially motivated by the original experiments of Fultz et al. (1959), where a heat source rotated around a cylinder with the radius R, except in our case thermal waves travel at the bottom and top and a mean temperature gradient was applied, as in Setup B of the previous part. In particular, the temperature distribution is linear in the radial r-direction and consists of one wave period in ϕ that travels counterclockwise: θ(ϕ, r, z = 0, t) = 0.5 r R cos(ϕ − 2πΩt) + 1 , θ(ϕ, r, z = H, t) = 0.5 r R cos(ϕ − 2πΩt) − 1 . Again, the mean temperature gradient, averaged over time, is the same as in classical RBC. The cell is shown in figure 7 (a). Furthermore, top and bottom plates are free-slip (∂u/∂n = 0) and no-slip conditions are applied at the sidewall (u = 0). All simulations are carried out for the parameters P r = 1 and the aspect ratio Γ ≡ R/H = 1. The rather large aspect ratio is a sacrifice, in return more simulations could be conducted and the parameter space in the region of interest is well resolved, as is shown in figure 7 (b).

Results
Previously, we have shown that travelling thermal waves generate a mean horizontal, or, synonymous, a zonal flow. The same can be observed in the cylindrical system, where a zonal flow now refers to non-vanishing azimuthal mean flow. In the following, we evaluate its strength and direction and discuss the results in the context of the 2D results. As no specific adjustments to the theoretical model have been done, from that point on, the model results are intended to serve mainly as references to the previous results. A brief remark beforehand: Evaluating the time and volume average of u ϕ proves problematic, as often flows are not purely pro-or retrograde. Therefore, rather than to give precise scaling laws, the primary purpose of the subsequent analysis is to explore the parameter space, demonstrate the overall strength of the zonal flows, find the most critical wave frequencies and determine the critical Ra above which the results deviate substantively from the predictions. Figure 8 shows (a) the total mean azimuthal momentum U ϕ V and (b) the value of U ϕ r,ϕ at the mid-height. As before, circles denote a retrograde, stars a prograde mean flow, and the solid lines are the 2D model solutions from Davey (1967), without modifications for no-slip walls. The obtained flows for small Ra 10 5 share distinct features with the 2D flows. The mean momentum converges to the asymptotic scalings, and, in fact, the data of figure 8 b collapse under a transformation with Ra remarkably good. For larger Ω, in particular Ω 10 −1 , the most flows are found to be directed prograde, even for Ra = 10 3 , which is different from the 2D case. And as in 2D, the flow structures reveal a transition in this Ω-region. As was discussed in section 3.1.2, the plane of the LSC drifts with the same speed as the TW (= Ω), if the TW speed is small compared to thermal diffusion speed, and the LSC breaks off from the TW at larger Ω, forming separate structures, acting on different timescales. It is in the regime of this Figure 9. For a fixed TW frequency Ω = 0.01. The azimuthally averaged mean azimuthal velocity Uϕ ϕ (top row) and the corresponding snapshots of the temperature θ (bottom row). As Ra increases, the core zonal flow becomes first stronger retrograde (Ra = 10 4 , 10 5 ), then switches its state to a prograde flow originated from the sidewall (Ra 10 6 ), while still increasing its strength (see colour bar). break-off above which a prograde flow is present. This process hints towards a similar mean flow instability, as discussed in §3.1.1, where the mean flow is now a slow LSC.
As Ra exceeds 10 5 , turbulent fluctuations increase and the data in figure 8 becomes increasingly scattered. The asymptotic scalings are hardly determinable, even though U x V ∼ Ω 1 for Ω → 0 appear still valid. The fluctuations can exceed their mean values, especially for small and large Ω. Despite the strong fluctuations, in regions of maximal zonal flow, i.e. Ω ≈ 10 −2 , the mean values are highly significant and can induce zonal flows of the same order of magnitude as the TW frequency, U ϕ V ≈ O(10 −2 ). Furthermore, similarly to the 2D case, in 3D, the zonal flows at high Ra are most of the time directed prograde, contrary to small Ra. From the vertical planes of the azimuthally and time averaged azimuthal velocity, shown in figure 9, the dominance of prograde motion at large Ra becomes more obvious. Moreover, these figures reveal a complex, inhomogeneous flow, with strong differential rotation and poloidal mean velocities.
In the following, we assess the contributing terms of the mean flow azimuthal momentum equation. For clarity, let us write the equation for u ϕ explicitly: (4.1) First, we consider how U ϕ changes in the vertical direction and, second, how it changes radially. Therefore, decomposing eq. (4.1) into its mean and fluctuation components, and averaging over ϕ and r gives the following balance: (4.2) Analysing the radial dependence, on the other side, averaging over ϕ and z gives P r Ra (4. 3) The rhs-terms of these equations are evaluated for Ω = 10 −2 , which are shown in figure  10. We ensured that in the simulations the data was averaged over an integer number of the TW periods, to prevent artifacts of the TW in the mean fields. The exact values are given in table 1. When we compare the individual mean velocities for (a) Ra = 10 3 and (b) Ra = 10 4 , it becomes clear that the mean field transport in both, vertical and radial, directions is rather negligible. Hence, the non-linear Reynolds stress sustains the mean zonal flow, just like in the 2D case for small Ra (see figure 4 a), and as expected (Stern 1959;Davey 1967). The small mean field contributions even reinforce the zonal flow, since the shape of the mean advection curves matches the shape of the Reynolds stress curve. Comparing further the vertical and the radial transport, we find that the former dominates the latter one by an order of magnitude. This proves that in this case the neglect of the radial currents, as suggested by Stern (1959), is justified, and therefore the mean momentum scalings (figure 8) match remarkably well with its 2D analogue (figure 3), and the difference in the prefactors can presumably be explained by the different velocity BCs.
The situation for larger Ra (figure 4 c-e) is vastly different. First, the problem becomes considerably three dimensional and the radial transport now reaches the same order of magnitude as the vertical transport (e.g. figure 10 c-e), which suggests that the validity of the 2D analogy at large Ra is no longer justified. Furthermore, the mean field advection contributions, which can be partially seen from figure 9, increase significantly. As a matter of fact, locally it can even exceed the Reynolds stress contributions. Furthermore, whereas for small Ra, vertical and radial momentum transports are present throughout the whole domain, at large Ra it becomes strongly confined to the boundaries. Especially the vertical transport peaks close to the top and bottom boundaries and is less pronounced in the center. The radial transport, on the other side, shows an interesting feature in the region 0.95 r R (figure 10 d,e). All terms are simultaneously positive, which causes an enhanced zonal transport close to the sidewall. This may explain why a prograde flow first appears close to the sidewall (figure 9, Ra = 10 6 ) and from there spreads further inwards (figure 9, Ra = 10 7 ).
Finally, we would like to illustrate the strength of the induced zonal flows on a concrete example. Assume an atmospheric boundary layer with a height ofĤ = 500m and a vertical temperature difference of ∆θ = 3 • . Given a mean temperature of 10 • , the material properties of air are approximately κ = 2.0 × 10 −5 m 2 /s, ν = 1.4 × 10 −5 m 2 /s and α = 3.6×10 −3 1/K. From that we find P r ≈ 0.7 and Ra ≈ 10 16 and the free-fall units u f f ≡ αgĤ∆θ ≈ 7m/s, t f f ≡Ĥ/u f f ≈ 70s. This system is exposed to a travelling thermal wave through the solar radiation with a period of 24h, or, in dimensionless units Vertical balance Ω ≈ 10 −3 . For simplicity, we say, the day and night difference is also about 3 • , which is likely to be a rather conservative estimate. Our study does not conclusively show, how the zonal flows scale up to Ra = 10 16 , but the results suggest a saturation at higher Ra, therefore we proceed using the maximum order of magnitude, which is U ϕ ≈ 10 −2 (for the given Ω it might be smaller). With these values, the thermal variation of the Earth's surface would induce a prevailing zonal flow of around 0.07m/s, or equivalently 0.3km/h. However, locally it could exceed this value (see figure 9) multiple times, therefore speeds of 1km/h are conceivable. Nevertheless, the variance of this estimate is rather high. Subsequent studies have to examine the influence of Ra, P r and the geometry, in order to make more confident statements about natural systems.

Conclusions
We have explored the original moving heat source problem by means of direct numerical simulations in 2D and 3D systems, for varying Rayleigh numbers Ra and travelling thermal wave (TW) frequency Ω. In the seminal works of Fultz et al. (1959) and Stern (1959), it was discovered that a system subjected to such a travelling wave generates Reynolds stresses, which induce a large scale mean horizontal, or equivalently zonal, flow directed counter to the propagating thermal wave. Therefore, in the first part, we revisited the theoretical model proposed by Davey (1967) and and found excellent agreement with the theory for low Ra flows, where even the absolute magnitude of the zonal flows is reproduced remarkably good. As Ra increases, the theoretical model overestimates the DNS data, which is consistent with the effects of higher order non-linear contributions (Whitehead 1972;Young et al. 1972;Hinch & Schubert 1971).
However, when an unstable mean temperature gradient is added to the system, the flows deviate substantially from the initial predictions and often reverse their direction to a prograde moving zonal flow. Such a behaviour was theorised before, as the result of a mean flow instability caused by the tilt of convection cells (Thompson 1970;Busse 1972Busse , 1983. Therefore, we have conducted a global linear stability analysis of a base flow near onset of convection and confirmed this hypothesis. The most unstable mode can give rise to a reverse of the horizontal velocity profile. Despite the strong plausibility, that this mean flow instability is the dominating mechanism at large Ra, the question remains open why prograde flow are more numerous than retrograde flows, while the mean flow instability suggests a spontaneous break of symmetry and therefore a more balanced distribution. In this context, it would be interesting to study in the future the interaction between the TW induced and convection rolls induced fields. In the second part we have examined the moving heat source problem in the context of a 3D cylindrical RBC system. The asymptotic scalings U ϕ V ∼ Ω 1 for Ω → 0 and U ϕ V ∼ Ω −4 for Ω → ∞ of the 2D theoretical model (Davey 1967) still hold in this system, especially at small Ra. An analysis of the vertical and radial momentum transport contributions suggests that the radial transport is negligible at small Ra, (which justifies a 2D approximation) but becomes relevant as Ra increases. Furthermore, again, large Ra are found predominantly inducing a prograde mean zonal flow. This gives more evidence that the prograde prevalence is likely not fully explained by the mean flow instability picture and further study's are required to explain its origin.
The studied problem is sufficiently general enough and can be extended to more complicated systems (Whitehead 1975;Shukla et al. 1981;Mamou et al. 1996). A more generalized theoretical framework already exist, which includes the influence of a basic stability and rotation (Stern 1971;Chawla & Purushothaman 1983), however, as this study showed, the theoretical models most often cannot fully explain the phenomena in convection dominated systems. Furthermore, the moving heat source problem might help to understand ubiquitous structures present in rotating systems. In rotating RBC systems, the flow structures near the sidewall (Favier & Knobloch 2020;Zhang et al. 2020) are similar to a certain extent to those structures accounting from the imposed TW.
Ultimately, this study also revealed that the estimate of the order of magnitudes is still afflicted with too large variances to make reliable statements about natural systems. A naive approach showed that atmospheric currents, caused by solar radiation and Earth's rotation, can actually generate prevailing zonal flows of about 1.0 km h . However, the variance of this estimate is rather high, it therefore is pivotal for subsequent studies to examine the sensitivities with Ra, P r and geometry in greater detail. setup, where the thermal wave travels only at the bottom, while the dimensionless top temperature equals zero.