1. 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 (Reference Halley1687) was probably the first to perceive that the periodic heating of the Earth's surface, due to the Earth's rotation, could be the reason for the occurrence of zonal winds in the atmosphere. Nearly three centuries later, experiments by Fultz et al. (Reference Fultz, Long, Owens, Bohan, Kaylor and Weil1959), 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 have appeared, which illuminated this phenomenon.
Thus, Stern (Reference Stern1959) repeated Fultz's experiments using a cylindrical annulus. His observations confirmed the previous result 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 (2-D) model, showing that the mean motion is maintained through the presence of the Reynolds stresses. Davey (Reference Davey1967) extended Stern's model and provided a 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, 
$\langle U_x\rangle _V$, with the characteristic frequency of the moving heat source 
$\varOmega$: 
$\langle U_x\rangle _V\sim \varOmega ^1$ for 
$\varOmega \rightarrow 0$ and 
$\langle U_x\rangle _V\sim \varOmega ^{-4}$ for 
$\varOmega \rightarrow \infty$. The topic gained further attention when Schubert & Whitehead (Reference Schubert and Whitehead1969) 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 (
$Pr$) 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 (
$Pr\ll 1$), by up to 4 times. This validated the linear analysis by Davey, who predicted the speed of the fluid to increase as 
$Pr$ 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 (Reference Whitehead1972), Young, Schubert & Torrance (Reference Young, Schubert and Torrance1972) and Hinch & Schubert (Reference Hinch and Schubert1971) 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 
$Pr$, as soon as convective processes come into play.
The preceding analysis certainly lacked the complexity of convective flows, and therefore Malkus (Reference Malkus1970), Davey (Reference Davey1967) and other authors anticipated that convective and shear instabilities could alter the entire character of the solution. In particular, Thompson (Reference Thompson1970) 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, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009) was studied theoretically by Busse (Reference Busse1972, Reference Busse1983) and Howard & Krishnamurti (Reference Howard and Krishnamurti1986), 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. Reference Goluskin, Johnston, Flierl and Spiegel2014), but only if the ratio of the horizontal to vertical extension of the domain is smaller than a certain value, see Wang et al. (Reference Wang, Verzicco, Lohse and Shishkina2020b) and Wang et al. (Reference Wang, Chong, Stevens and Lohse2020a). Also, several studies examined the effects of time-dependent sinusoidal perturbations in RBC. Venezian (Reference Venezian1969) showed that the onset of convection can be advanced or delayed by modulation, while Yang et al. (Reference Yang, Chong, Wang, Verzicco, Shishkina and Lohse2020) and Niemela & Sreenivasan (Reference Niemela and Sreenivasan2008) 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 set-up with a 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 TW problem is merely of academic interest or, indeed, of practical relevance in the generation of geo- and astrophysical zonal flows (Yano, Talagrand & Drossart Reference Yano, Talagrand and Drossart2003; Maximenko, Bang & Sasaki Reference Maximenko, Bang and Sasaki2005; Nadiga Reference Nadiga2006). For this purpose, in § 4, we complement our analysis with thorough 3-D DNS. For the sake of generality, we choose a classical RBC set-up. Ultimately, we analyse the absolute angular momentum in 3-D flows (respectively, horizontal velocity in 2-D flows) and provide insight into the mean flow structures.
2. Methods
2.1. Direct numerical simulations
The governing equations in the Oberbeck–Boussinesq approximation for the dimensionless velocity 
$\boldsymbol {u}$, temperature 
$\theta$ and pressure 
$p$ read as follows:
\begin{gather} {\partial }\boldsymbol{u}/{\partial } t+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}+\boldsymbol{\nabla} {p}= \sqrt{Pr/Ra} \nabla^2 \boldsymbol{u}+ {\theta}\boldsymbol{e}_z , \end{gather}
\begin{gather}{\partial }{\theta}/{\partial } t+\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} {\theta}=1/\sqrt{Pr Ra} \nabla^2 {\theta}, \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} =0. \end{gather}
Here, 
$t$ denotes time and 
$\boldsymbol {e}_z$ the unit vector in the vertical direction. The equations have been non-dimensionalised using the free-fall velocity 
$u_{ff}\equiv (\alpha g \Delta \hat {H})^{1/2}$, the free-fall time 
$t_{ff}\equiv \hat {H}/u_{ff}$, 
$\varDelta$ the amplitude of the thermal TW and 
$\hat {H}$ the cell height. The dimensionless parameters 
$Ra$, 
$Pr$ and the aspect ratio 
$\varGamma$ are defined by
\begin{equation} Ra\equiv \alpha g \Delta \hat{H}^3/(\kappa\nu), \quad Pr\equiv\nu/\kappa, \quad \varGamma\equiv \hat{L}/\hat{H}, \end{equation}
where 
$\hat {L}$ is the length of the domain, 
$\nu$ is the kinematic viscosity, 
$\alpha$ the isobaric thermal expansion coefficient, 
$\kappa$ the thermal diffusivity and 
$g$ the acceleration due to gravity. This set of equations is solved using the finite-volume code goldfish (Shishkina et al. Reference Shishkina, Horn, Wagner and Ching2015; Kooij et al. Reference Kooij, Botchev, Frederix, Geurts, Horn, Lohse, van der Poel, Shishkina, Stevens and Verzicco2018), which employs a fourth-order discretisation scheme in space and a third-order Runge–Kutta, or, alternatively, an Euler-leapfrog scheme in time. The code runs on rectangular and cylindrical domains and has been advanced for a 2-D pencil decomposition for highly parallel usage. The spatial grid resolution of the simulations was chosen according to the minimum resolution requirements of Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010). A stationary state is ensured by monitoring the volume-averaged, the wall-averaged and the kinetic dissipation based Nusselt numbers.
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^\prime$, decomposing 
$u$ into its mean part 
$U$ and fluctuating part 
$u^\prime$. Unless specifically stated, time averages are carried out over a long period of time, however, in § 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 
$\langle \cdot \rangle$, followed by the respective direction of the average, e.g. 
$\langle \cdot \rangle _x$ denotes an average in 
$x$; 
$\langle \cdot \rangle _V$ denotes a volume average. And ultimately, the velocity vector definitions 
$\boldsymbol {u} \equiv (u_x,u_y,u_z) \equiv (u,v,w)$ are used interchangeably.
2.2. Theoretical model
Already the earliest models proposed by Stern (Reference Stern1959) and Davey (Reference Davey1967) gave a considerable good understanding of the moving heat source problem. Although there are more complex models (Stern Reference Stern1971) based on adding higher-order nonlinear contributions (Hinch & Schubert Reference Hinch and Schubert1971; Busse Reference Busse1972; Whitehead Reference Whitehead1972; Young et al. Reference Young, Schubert and Torrance1972), 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 linearised 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:
\begin{equation} \sqrt{Pr/Ra} {\partial_z^2} \langle U \rangle_x = \partial_z \langle \overline{u^\prime w^\prime} \rangle_x + \langle W \partial_z U \rangle_x. \end{equation}
Evidently, a mean zonal flow 
$\langle U \rangle _x$ is maintained by the momentum transport due to the Reynolds stress component 
$\overline {u^\prime w^\prime }$ and by mean advection through 
$W \partial _z U$. The theory further advances by assuming that no vertical mean flow exists (
$W=0$), which reduces (2.4) to the balance between viscous mean diffusion and Reynolds stress diffusion. Furthermore, by neglecting convection and variations in 
$x$, the linearised equations can be written as a set of ordinary differential equations, that can be solved sequentially to find 
$u^\prime$ and 
$w^\prime$ and ultimately the Reynolds stress term 
$\overline {u^\prime w^\prime }$. This procedure is shown in appendix A. Given the Reynolds stress field, (2.4) has to be integrated twice to obtain the mean zonal flow 
$U(z)$. Integrating that profile again finally gives the total mean zonal flow 
$\langle U \rangle _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 (Reference Davey1967), the limiting relations can be calculated explicitly
\begin{gather} \langle U_x \rangle_V ={-}\frac{\rm \pi}{2}\frac{ k^3 Ra^2 Pr^{{-}2}(Pr+1)}{ 12!} \varOmega + {O}(\varOmega^3) \quad \text{for}\ \varOmega \rightarrow 0, \end{gather}
\begin{gather}\langle U_x \rangle_V ={-}\frac{k^3 Ra^{{-}1/2} Pr^{{-}3/2}}{ 256 {\rm \pi}^4 (Pr+1)} \varOmega^{{-}4} + {O}(\varOmega^{{-}9/2}) \quad \text{for}\ \varOmega \rightarrow \infty, \end{gather}
where the horizontal TW, 
$\theta (x,t)=0.5 \cos (kx-2{\rm \pi} \varOmega t)$, is applied to the bottom and top plate. 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 \partial _z U$-contribution is no longer justified.
3. Two-dimensional convective system
As described by Stern (Reference Stern1959), the generation of a laminar zonal flow by a TW can be successfully explained in a 2-D system, which makes it a good starting point. The temperature boundary conditions (BCs) are time and space dependent,
\begin{gather} \theta(x,z=0,t) = 0.5\left[\cos(x-2{\rm \pi}\varOmega t) + \Delta \theta \right], \end{gather}
\begin{gather}\theta(x,z=H,t) = 0.5\left[\cos(x-2{\rm \pi}\varOmega t) - \Delta \theta \right]. \end{gather}
Here, 
$\varOmega$ indicates the temporal frequency of the TW in free-fall time units. For example, 
$\varOmega =10^{-1}$ describes a wave with a period of 
$10$ free-fall time units 
$\tau _{ff}$, and 
$\Delta \theta$ is introduced as a control parameter for the strength of the mean temperature gradient.
In the following, two different set-ups are considered. In set-up A (figure 1a) – the one originally examined by Davey (Reference Davey1967) – no mean temperature gradient exists (
$\Delta \theta =0$) and the top and bottom plate temperatures are equal, whereas in set-up B (figure 1b) a mean, unstable temperature gradient is applied (
$\Delta \theta =1$). For simplicity, the mean temperature gradient is set equal to the amplitude of the thermal wave. In this set-up, effects of convection are expected to become dominant. Averaged over time, this set-up 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 length 
$L=2{\rm \pi}$ and height 
$H=1$. In upcoming studies, one might introduce a second Rayleigh number based on the mean temperature gradient (as in RBC), namely 
$Ra_{\Delta \theta } \equiv \alpha g \Delta \theta \hat {H}^3/(\kappa \nu )$. However, in this work the connection to 
$Ra$ is simply 
$Ra_{\Delta \theta }=0$ for set-up A and 
$Ra_{\Delta \theta }=Ra$ for set-up B.
Figure 1. Sketch of the 2-D numerical set-up. The colour represents the dimensionless temperature distribution for the purely conductive cases (
$\varOmega =0.1$). The thermal wave is imposed at the top and bottom plates, propagating to the right, in the positive 
$x$-direction. (a) Set-up A: no mean temperature gradient is imposed between the top and the bottom. (b) Set-up B: with (unstably stratified) mean temperature gradient, as in RBC. The figure shows temperature snapshots, while the time-averaged conduction temperature field depends linearly on 
$z$.
The overall focus in this study lies on variations of the zonal flow with 
$Ra$ and 
$\varOmega$. Thus, the parameter space spans 
$10^3 \leq Ra \leq 10^7$ and 
$10^{-4}\leq \varOmega \leq 10^0$, while the aspect ratio and Prandtl number are kept constant (
$\varGamma =2{\rm \pi}$, 
$Pr=1$). Exemplary temperature fields at a fixed 
$\varOmega =0.1$ are shown in figure 2.
Figure 2. Snapshots of the temperature field 
$\theta$ at a fixed TW speed 
$\varOmega =0.1$ (propagating to the right). (a) Set-up A and (b) set-up B. The plumes in set-up B travel either retrograde or prograde (see supplementary movies available at https://doi.org/10.1017/jfm.2020.1186).
3.1. Results
The theoretical model, as presented in appendix A, aims to explain the generation of the total mean momentum 
$\langle U_x \rangle _V$ for a given 
$Ra$ and wave frequency 
$\varOmega$. 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 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 time scale rather than the free-fall time scale 
$\varOmega _{\kappa ,max} = \varOmega /\sqrt {RaPr} \approx 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.
Figure 3. Mean velocity of the zonal flow vs. the wave frequency 
$\varOmega$ for 
$Ra = 10^3$ (blue), 
$10^4$ (orange), 
$10^5$ (green), 
$10^6$ (red) and 
$10^7$ (black). 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 (Reference Davey1967). (a) Set-up A and (b) set-up B.
We begin our discussion with the results of set-up A, shown in figure 3(a). The theoretical model by Davey (Reference Davey1967), 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 (Reference Whitehead1972), Young et al. (Reference Young, Schubert and Torrance1972) and Hinch & Schubert (Reference Hinch and Schubert1971) who observed that corrections of higher-order nonlinear 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 
$\varOmega$ are well predicted by the model, the large 
$Ra$ flows are mostly over predicted (except 
$Ra=10^7$ and 
$\varOmega =0.1$, the only flow of that set-up that becomes truly turbulent, despite similar initial conditions). Presumably, even more important is that some of the flows for 
$Ra \geq 10^5$ exhibit a positive/prograde mean flow, indicated by a star symbol, which is especially prevalent at small 
$\varOmega$.
Turning the focus to set-up 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\approx 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 \ge 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-hand side of (2.4), which are presented in figure 4. The model neglects mean advection, it only captures contributions of 
$\overline {u^\prime w^\prime }$. 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 \partial _z U$ starts to take over. The shape of the mean flow advection curve is antiphase to the Reynolds stress curve and contributes the most. This explains the reversal of the mean flow, from retrograde in figure 4(a) to prograde in figure 4(b).
Figure 4. Mean profiles of Reynolds stress vs. mean flow advection contribution for 
$Ra=10^4$ and 
$\varOmega =0.1$ for (a) 
${\Delta \theta =0}$ and (b) 
${\Delta \theta =1}$, where mean advection dominates. The flow in (a) moves retrograde due to the Reynolds stress contribution, while the flow in (b) shows a prograde mean flow (
$\sqrt {Pr/Ra} {\partial _z^2} \langle U \rangle _x = \partial _z \langle \overline {u^\prime w^\prime } \rangle _x + \langle W \partial _z U \rangle _x$).
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 free-fall time rather than with the diffusive time and the maximum is found in the region 
$0.01\leq \varOmega \leq 0.1$. This is important, since natural flows often fall within this parameter range. We show this in the context of the Earth's atmosphere in § 4. Finally the instantaneous fields most often show three plumes (figure 2b), 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 pre-existing 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.
In figure 5 we show the total kinetic 
$E_{tot}=(\langle \overline {u_x^2+u_z^2} \rangle _V^{1/2})/2$ and horizontal (zonal flow) kinetic energy 
$E_{hor}=(\langle \overline {u_x^2}\rangle _V^{1/2})/2$ in order to elucidate the energetic impact of the present zonal flows and to evaluate the strength of the vertical and horizontal motions. Set-up A (a) is clearly dominated by the horizontal kinetic energy throughout the whole parameter range. For 
$\varOmega >10^{-1}$, the kinetic energy drops close to zero and the temperature is transported by conduction only above this limit. However, before the kinetic energy drops, the curves show an energy enhancement. The location of the energy maximum coincides with the maximum of the zonal flow (figure 3a), which indicates that the zonal flows can have a significant imprint in the energy of the system. Likewise, set-up B (figure 5b) is also dominated by horizontal kinetic energy. However, obviously for larger 
$Ra$ and larger 
$\varOmega$, the magnitude of the vertical kinetic energy becomes increasingly important. This further supports that the neglect of the vertical velocity component 
$W$ is eventually no longer justified for these parameter regimes.
Figure 5. Total kinetic energy 
$E_{tot}$ (black) and horizontal kinetic energy 
$E_{hor}$ (blue) for (a) set-up A and (b) set-up B.
3.1.1. 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. Set-up B and 
$\varOmega =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 3b). Thus, a simulation was initiated at 
$Ra=1000$ and then 
$Ra$ was progressively increased by 
$1000$ each time after a steady state had settled. The time evolution of the total mean zonal flow is given in figure 6(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 6(b), flattens. Ultimately, at 
$Ra\geq 4000$ this profile flips over and the total zonal flow turns into a prograde state.
Figure 6. Path from a retrograde flow to a prograde flow. (a) Time evolution of the mean zonal flow; 
$Ra$ was increased stepwise. For 
$Ra\geq 3000$ convection rolls form; for 
$Ra\geq 4000$ the rolls tilt significantly, and the mean zonal flow becomes positive. (b) The 
$u_x$ profiles for 
$Ra=1000, 3000, 5000$. (c) Mean flow extracted at 
$Ra=3000$ (averaged over one TW period). (d) Result of the global stability analysis for the mean flow of (c), that becomes unstable for 
$Ra\geq 4000$ to tilted convection rolls.
As we have shown in the preceding analysis (figure 4b), 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\partial _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 
$\partial _z U$. A possible mechanism, even discussed in the context of the moving heat source problem, was described by Thompson (Reference Thompson1970) and theoretically analysed by Busse (Reference Busse1972), 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 (Reference Busse1972) 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 6(a) at 
$Ra=3000$ coincides with the observed onset of convection, which is slightly delayed compared to classical, unmodulated RBC (
$Ra_c\approx 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 time scales, which results in the base flow, as shown in figure 6(c). Based on this base flow, a linear, temporal stability analysis of the full 2-D linearised Navier–Stokes equations was conducted. Details therefore are given in 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 6(d). The growth rate of it is 
$\sigma \approx 0+0.2{\rm i}$, 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 (Reference Thompson1970), in the sense that the mode induces a mean shear flow (see profile in figure 6d). However, rather than the ‘pure’ shear flows as presented by Thompson (Reference Thompson1970) and Busse (Reference Busse1972, Reference Busse1983) with a vanishing total net momentum when integrated vertically, the fluctuation profile found in our study (figure 6d on the right) shows a more directed flow, negative in the vicinity of the plates and stronger positive in the centre. 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 6(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 (Reference Krishnamurti and Howard1981), it was already remarked that the mean flow transition is caused by a spontaneous symmetry breaking 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 6(d). And finally, even though in figure 6(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.
3.1.2. 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 2-D system and, already ahead of the next part, in the 3-D cylindrical system. In addition, movies are provided as supplementary material.
In general, in two dimensions, as can be seen from figure 2, the temperature field is either symmetric around the horizontal mid-plane (set-up A), or not; in this case there exist plumes (set-up B). In the latter case, there are usually three up- and three down-welling plumes identifiable. In the 3-D case, the flow consists of rising and falling plumes, which together form a large scale circulation (LSC). If the TW propagates slowly (small 
$\varOmega$), the plumes (two dimensions) or respectively the LSC plane (three dimensions) drift with the same speed as the TW and both structures appear to be connected. However, as 
$\varOmega$ increases and, hypothetically, the TW time scale 
$\tau _\varOmega$ becomes small compared to thermal diffusion 
$\tau _\kappa$ (
$\tau _\kappa /\tau _\varOmega = \sqrt {Pr Ra}\varOmega$), the plumes (two dimensions) or LSC (three dimensions) ‘break-off’ from the TW, forming two separate structures, acting on different time scales.
Figure 7 shows the space–time structures of the temperature field, evaluated at mid-height, and in the 3-D case at mid-height and near the sidewall. The structures at mid-height 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 
$\varOmega$ parameters, (ii) is found for large 
$Ra$ and large 
$\varOmega$, if no mean temperature is present and (iii) exists in strongly convection dominated flows for large 
$Ra$ and large 
$\varOmega$, especially if a mean temperature gradient is present. Furthermore, it is striking that temperatures between the left and right regions in the vicinity of the plumes centre (hottest or coldest regions in figure 7) do not necessarily fill with the same temperature (c). This gives further evidence of a mean flow instability, as it features similarities of the temperature field of the unstable mode given in figure 6(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, nonlinear saturation.
Figure 7. Evolution of the temperature at mid-height 
$z=H/2$, for (a–c) 2-D flows and (d,e) 3-D RBC (
$r=0.99R$) flows; (a) 
$\varOmega =0.01$, set-up A, (b) 
$\varOmega =0.1$, set-up A, (c) 
$\varOmega =0.1$, set-up B, (d) 
$\varOmega =0.01$, 3-D RBC and (e) 
$\varOmega =0.1$ 3-D RBC. All panels show 
$Ra=10^7$. The black solid line in each plot indicates the TW speed.
4. Three-dimensional convective systems
The preceding part, as most of the existing literature, is confined to 2-D flows. Now we will discuss the moving heat source problem in the context of more complicated 3-D convective flows. In general, TW solutions are common amongst 3-D convective systems. Bensimon et al. (Reference Bensimon, Croquette, Kolodner, Williams and Surko1990), Kolodner, Bensimon & Surko (Reference Kolodner, Bensimon and Surko1988) and Kolodner & Surko (Reference Kolodner and Surko1988) observed convection rolls propagating azimuthally in a large aspect ratio annulus near the onset of convection. Their drift velocity was 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 TW solution in RBC systems are the spiral patterns found in large aspect ratio cells (Bodenschatz et al. Reference Bodenschatz, de Bruyn, Ahlers and Cannell1991; Bodenschatz, Pesch & Ahlers Reference Bodenschatz, Pesch and Ahlers2000). These spirals are rotating in either direction, although corotating spirals are more numerous (Cross & Tu Reference Cross and Tu1995), and are known to be coupled with an azimuthal mean flow (Decker, Pesch & Weber Reference Decker, Pesch and Weber1994). Furthermore, in rotating systems, TW structures are quite common (Knobloch & Silber Reference Knobloch and Silber1990). These structures are strongly geometry dependent (Wang et al. Reference Wang, Ma, Chen and Sun2012) and known to induce mean zonal flows that propagate pro- and retrograde (Zhang et al. Reference Zhang, van Gils, Horn, Wedi, Zwirner, Ahlers, Ecke, Weiss, Bodenschatz and Shishkina2020).
Despite the vast literature on these phenomena, quantitative data on mean flows that are induced by external travelling thermal waves in 3-D flows seem to be rare. Therefore our main goal in this part is to gain insight into 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 DNS.
4.1. Numerical set-up: cylindrical RBC (
$Pr=1$)
The set-up is essentially motivated by the original experiments of Fultz et al. (Reference Fultz, Long, Owens, Bohan, Kaylor and Weil1959), where a heat source rotated around a cylinder with the radius 
$R$ (diameter 
$D$), except, in our case, thermal waves travel at the bottom and top and a mean temperature gradient was applied, as in set-up B of the previous part. In particular, the temperature distribution is linear in the radial 
$r$-direction and consists of one wave period in 
$\varphi$ that travels counterclockwise
\begin{gather} \theta(\varphi,r,z=0,t) = 0.5\left[\frac{r}{R} \cos(\varphi - 2{\rm \pi} \varOmega t) + 1\right], \end{gather}
\begin{gather}\theta(\varphi,r,z=H,t) = 0.5\left[\frac{r}{R} \cos(\varphi - 2{\rm \pi} \varOmega t) - 1\right]. \end{gather}
Again, the mean temperature gradient, averaged over time, is the same as in classical RBC. The cell is shown in figure 8(a). Furthermore, top and bottom plates are free slip (
$\partial \boldsymbol {u}/ \partial \boldsymbol {n} = 0$) and no-slip conditions are applied at the sidewall (
$\boldsymbol {u}=0$). All simulations are carried out for the parameters 
$Pr=1$ and the aspect ratio 
$\varGamma \equiv D/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 shown in figure 8(b).
Figure 8. (a) Sketch of the cylindrical domain and imposed TW. (b) Studied parameter space. The mesh sizes 
$n_r \times n_\varphi \times n_z$ of the DNS are 
$48 \times 130 \times 98$ for 
$Ra=10^3$, 
$96 \times 260 \times 196$ for 
$Ra=10^4,10^5,10^6$ and 
$128 \times 342 \times 256$ for 
$Ra=10^7$.
4.2. Results
Previously, we have shown that travelling thermal waves generate a mean horizontal, or, synonymously, 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 2-D results. As no specific adjustments to the theoretical model have been made, from this 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_\varphi$ proves problematic, as often flows are not purely pro- or retrograde. Therefore, rather than 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 and find the most critical wave frequencies and determine the critical 
$Ra$ above which the results deviate substantively from the predictions.
Figure 9(a) shows the total mean azimuthal momentum 
$\langle U_\varphi \rangle _V$ and figure 9(b) shows the value of 
$\langle U_\varphi \rangle _{r,\varphi }$ at the mid-height. As before, circles denote a retrograde, stars a prograde mean flow and the solid lines are the 2-D model solutions from Davey (Reference Davey1967), without modifications for no-slip walls. The obtained flows for small 
$Ra\leq 10^5$ share distinct features with the 2-D flows. The mean momentum converges to the asymptotic scalings, and, in fact, the data of figure 9(b) collapse under a transformation with 
$Ra$ remarkably well. For larger 
$\varOmega$, in particular 
$\varOmega \geq 10^{-1}$, the most flows are found to be directed prograde, even for 
$Ra=10^3$, which is different from the 2-D case. And as in two dimensions, the flow structures reveal a transition in this 
$\varOmega$-region. As was discussed in § 3.1.2, the plane of the LSC drifts with the same speed as the TW (
$=\varOmega$), if the TW speed is small compared to thermal diffusion speed, and the LSC breaks off from the TW at larger 
$\varOmega$, forming separate structures, acting on different time scales. It is in the regime of this 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.
Figure 9. (a) Time- and volume-averaged zonal flow as a function of the heat source frequency 
$\varOmega$, (b) zonal flow at mid-height for 3-D RBC data; 
$Ra =10^3$ (blue), 
$10^4$ (orange), 
$10^5$ (green), 
$10^6$ (red) and 
$10^7$ (black). 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 (Reference Davey1967).
As 
$Ra$ exceeds 
$10^5$, turbulent fluctuations increase and the data in figure 9 become increasingly scattered. The asymptotic scalings are hardly determinable, even though 
$\langle U_x\rangle _V\sim \varOmega ^1$ for 
$\varOmega \rightarrow 0$ appears still valid. The fluctuations can exceed their mean values, especially for small and large 
$\varOmega$. Despite the strong fluctuations, in regions of maximal zonal flow, i.e. 
$\varOmega \approx 10^{-2}$, the mean values are highly significant and can induce zonal flows of the same order of magnitude as the TW frequency, 
$\langle U_\varphi \rangle _V \approx 10^{-2}$. Furthermore, similarly to the 2-D case, in three dimensions, 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 10, 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.
Figure 10. For a fixed TW frequency 
$\varOmega = 0.01$. The azimuthally averaged mean azimuthal velocity 
$\langle U_\varphi \rangle _\varphi$ (a–e) and the corresponding snapshots of the temperature 
$\theta$ (f–j). As 
$Ra$ increases, the core zonal flow becomes first stronger retrograde (
$Ra=10^4,10^5$), then switches its state to a prograde flow originating from the sidewall (
$Ra \geq 10^6$), while still increasing its strength (see colour bar).
4.2.1. Vertical and radial momentum transport
In the following, we assess the contributing terms of the mean flow azimuthal momentum equation. For clarity, let us write the equation for 
$u_\varphi$ explicitly
\begin{align} &\partial_t u_\varphi + \frac{1}{r} \frac{\partial r u_\varphi u_r}{\partial r} + \frac{1}{r} \frac{\partial u_\varphi u_\varphi}{\partial \varphi} + \frac{\partial u_\varphi u_z}{\partial z}\nonumber\\ &\quad =-\frac{1}{r}\frac{\partial p}{\partial \varphi} + \sqrt{\frac{Pr}{Ra}} \left[ \frac{1}{r} \frac{\partial}{\partial r}\left( r \frac{\partial u_\varphi}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 u_\varphi}{\partial \varphi^2} +\frac{\partial^2 u_\varphi}{\partial z^2} - \frac{u_\varphi}{r^2} + \frac{2}{r^2}\frac{\partial u_r}{\partial \varphi}\right]. \end{align} First, we consider how 
$U_\varphi$ changes in the vertical direction and, second, how it changes radially. Therefore, decomposing (4.3) into its mean and fluctuating components, and averaging over 
$\varphi$ and 
$r$ gives the following balance:
\begin{align} &\sqrt{\frac{Pr}{Ra}} \left( {\frac{\partial^2 \langle U_\varphi \rangle_{r,\varphi}}{\partial z^2}} -\frac{\langle U_\varphi \rangle_{r,\varphi}}{r^2} \right)\nonumber\\ &\quad = \frac{\partial \langle \overline{u_\varphi^\prime u_z^\prime} \rangle_{r,\varphi}}{\partial z} + \frac{\partial \langle U_\varphi U_z \rangle_{r,\varphi}}{\partial z} + \left\langle \frac{ \overline{u_\varphi^\prime u_r^\prime }}{r}\right\rangle_{r,\varphi} + \left\langle \frac{U_r U_\varphi }{r}\right\rangle_{r,\varphi}. \end{align} Analysing the radial dependence, on the other hand, averaging over 
$\varphi$ and 
$z$ gives
\begin{align} &\sqrt{\frac{Pr}{Ra}} \left( \frac{1}{r^2} {\frac{\partial^2 \langle U_\varphi \rangle_{\varphi,z}}{\partial r^2}} -\frac{\langle U_\varphi \rangle_{\varphi,z}}{r^2} \right)\nonumber\\ &\quad = \frac{1}{r}\frac{\partial r\langle \overline{u_\varphi^\prime u_r^\prime} \rangle_{\varphi,z}}{\partial r} + \frac{1}{r}\frac{\partial r\langle U_\varphi U_r \rangle_{\varphi,z}}{\partial r} + \left\langle \frac{ \overline{u_\varphi^\prime u_r^\prime }}{r}\right\rangle_{\varphi,z} + \left\langle \frac{U_r U_\varphi }{r}\right\rangle_{\varphi,z}. \end{align} The right-hand side terms of these equations are evaluated for 
$\varOmega =10^{-2}$, which are shown in figure 11. We ensured that, in the simulations, the data were averaged over an integer number of TW periods, to prevent artefacts of the TW in the mean fields (the exact time values can be found in the supplementary material). 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 nonlinear Reynolds stress sustains the mean zonal flow, just as in the 2-D case for small 
$Ra$ (see figure 4a), as expected (Stern Reference Stern1959; Davey Reference Davey1967). 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 radial transports, 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 (Reference Stern1959), is justified, and therefore the mean momentum scalings (figure 9) match remarkably well with their 2-D analogue (figure 3), and the difference in the prefactors can presumably be explained by the different velocity BCs.
The situation for larger 
$Ra$ (figure 11c–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 11c–e), which suggests that the validity of the 2-D analogy at large 
$Ra$ is no longer justified. Furthermore, the mean field advection contributions, which can be partially seen from figure 10, increase significantly. As a matter of fact, locally, the mean field advection 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. In particular, the vertical transport peaks close to the top and bottom boundaries and is less pronounced in the centre. The radial transport, on the other hand, shows an interesting feature in the region 
$0.95 \leq r/R \leq 1$ (figure 11d,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 10, 
$Ra = 10^6$) and, from there, spreads further inwards (figure 10, 
$Ra = 10^7$).
4.2.2. Sensitivity to the BCs and aspect ratio
The systems studied in this paper allow many variations of the velocity and temperature boundary conditions as well as geometrical characteristics of the system. Discussing all of them goes beyond the scope of a single study. Nevertheless, in order to provide some preliminary intuition, we examine selected variations and their effects on the generation of the zonal flows. We do this for a single baseline simulation at 
$Ra=10^5$ and 
$\varOmega =10^{-1}$. The mean angular momentum profiles are shown in figure 12.
Figure 12. Mean angular momentum profile for 
$Ra=10^5$ and 
$\varOmega =10^{-1}$. The curves show the effects of different imposed BCs and aspect ratios: baseline simulation (black, free-slip BCs, 
$\theta \sim r/R$ and 
$\varGamma =1$), sinusoidal radial temperature BCs (blue, free-slip BCs, 
$\theta \sim sin({\rm \pi} r/R)$ and 
$\varGamma =1$), no-slip (red, no-slip BCs, 
$\theta \sim r/R$ and 
$\varGamma =1$), 
$\varGamma =0.2$ (yellow, free-slip BCs, 
$\theta \sim r/R$ and 
$\varGamma =0.2$) and 
$\varGamma =2$ (green, free-slip BCs, 
$\theta \sim r/R$ and 
$\varGamma =2$).
First, we consider the effects of the aspect ratio. From classical RBC, it is known that zonal flow properties depend strongly on 
$\varGamma$ (Wang et al. Reference Wang, Chong, Stevens and Lohse2020a). In our case, a decrease of the aspect ratio from 
$\varGamma =1$ to 
$\varGamma =0.2$ (slender cell) weakens the zonal flow considerably by a factor of 
$100$. Furthermore, the zonal flow becomes confined to the top and bottom plates, while no zonal flow is observed in the centre of the cell. On the other hand, increasing the aspect ratio to 
$\varGamma =2$ has only minor impact on the zonal flow. We must note that for the case of 
$\varGamma =0.2$, convection has yet not started and subsequent studies would be necessary to conclusively elucidate on the aspect ratio dependence.
The effects of the BC variations on the formation of zonal flows can be formulated as follows. No-slip conditions at the top and bottom plates lead to a slightly weaker, but qualitatively similar zonal flow. Likewise, replacing the linear radial temperature distribution at the plates by a sinusoidal distribution (
$\theta \sim \sin ({\rm \pi} r/R)$) shows still a qualitatively similar angular momentum profile, although the strength of the zonal flow in the centre of the cell increases by a factor of approximately 
$1.5$. This indicates that the system is rather sensitive to variations of the temperature BCs.
4.3. Example: atmospheric boundary layer
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 
$\hat {H}=500\ \text{m}$ and a vertical temperature difference of 
$\Delta T=3\,^\circ \textrm {C}$. Given a mean temperature of 
$10\,^\circ \textrm {C}$, the material properties of air are approximately 
$\kappa =2.0 \times 10^{-5}\ \textrm {m}^2\,\textrm {s}^{-1}$, 
$\nu =1.4 \times 10^{-5}\ \textrm {m}^2\,\textrm {s}^{-1}$ and 
$\alpha =3.6 \times 10^{-3}\ \textrm {K}^{-1}$. From that, we find 
$Pr\approx 0.7$ and 
$Ra\approx 10^{16}$ and the free-fall units 
$u_{ff} \equiv \sqrt {\alpha g \hat {H} \Delta \theta }\approx 7\ \textrm {m}\,\textrm {s}^{-1}$, 
$t_{ff} \equiv \hat {H}/u_{ff} \approx 70\ \textrm {s}$. This system is exposed to a travelling thermal wave through the solar radiation with a period of 
$24$ h, or, in dimensionless units 
$\varOmega \approx 10^{-3}$. For simplicity, we say, the day and night difference is also approximately 
$3\,^\circ \textrm {C}$, 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_\varphi \approx 10^{-2}$ (for the given 
$\varOmega$ it might be smaller). With these values, the thermal variation of the Earth's surface would induce a prevailing zonal flow of around 
$0.07\ \textrm {m}\,\textrm {s}^{-1}$, or equivalently 
$0.3\ \textrm {km}\,\textrm {h}^{-1}$. However, locally, it could exceed this value (see figure 10) multiple times, therefore speeds of 
$1\ \textrm {km}\,\textrm {h}^{-1}$ are conceivable. Nevertheless, the variance of this estimate is rather high. Subsequent studies have to examine the influence of 
$Ra,Pr$ and the geometry, in order to make more confident statements about natural systems.
5. Conclusions
We have explored the original moving heat source problem by means of DNS in 2-D and 3-D systems, for varying Rayleigh numbers 
$Ra$ and travelling thermal wave frequency 
$\varOmega$. In the seminal works of Fultz et al. (Reference Fultz, Long, Owens, Bohan, Kaylor and Weil1959) and Stern (Reference Stern1959), it was discovered that a system subjected to such a TW 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 (Reference Davey1967) and found excellent agreement with the theory for low 
$Ra$ flows, where even the absolute magnitude of the zonal flows is reproduced remarkably well. As 
$Ra$ increases, the theoretical model overestimates the DNS data, which is consistent with the effects of higher-order nonlinear contributions (Hinch & Schubert Reference Hinch and Schubert1971; Whitehead Reference Whitehead1972; Young et al. Reference Young, Schubert and Torrance1972).
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 to be the result of a mean flow instability caused by the tilt of convection cells (Thompson Reference Thompson1970; Busse Reference Busse1972, Reference Busse1983). 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 as to 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 3-D cylindrical RBC system. The asymptotic scalings 
$\langle U_\varphi \rangle _V\sim \varOmega ^1$ for 
$\varOmega \rightarrow 0$ and 
$\langle U_\varphi \rangle _V\sim \varOmega ^{-4}$ for 
$\varOmega \rightarrow \infty$ of the 2-D theoretical model (Davey Reference Davey1967) 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 2-D approximation), but becomes relevant as 
$Ra$ increases. Furthermore, again, large 
$Ra$ is found to predominantly induce 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 studies are required to explain its origin.
The studied problem is sufficiently general and can be extended to more complicated systems (Whitehead Reference Whitehead1975; Shukla et al. Reference Shukla, Yu, Rahman and Spatschek1981; Mamou et al. Reference Mamou, Robillard, Bilgen and Vasseur1996). A more generalised theoretical framework already exists, which includes the influence of a basic stability and rotation (Stern Reference Stern1971; Chawla & Purushothaman Reference Chawla and Purushothaman1983), 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 the ubiquitous structures present in rotating systems. In rotating RBC systems, the flow structures near the sidewall (Favier & Knobloch Reference Favier and Knobloch2020; Zhang et al. Reference Zhang, van Gils, Horn, Wedi, Zwirner, Ahlers, Ecke, Weiss, Bodenschatz and Shishkina2020) are similar to a certain extent to those structures due to the imposed TW.
Ultimately, this study also revealed that the estimates of the order of magnitudes are 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 the Earth's rotation, can actually generate prevailing zonal flows of approximately 
$1.0\ \textrm {km}\,\textrm {h}^{-1}$. However, the variance of this estimate is rather high, it therefore is pivotal for subsequent studies to examine the sensitivities with 
$Ra,Pr$ and the geometry in greater detail.
Supplementary material and movies
Supplementary material and movies are available at https://doi.org/10.1017/jfm.2020.1186.
Acknowledgements
The authors would like to thank the Max-Planck HPC Teams in Göttingen and Munich for their generous technical support and additional computational resources.
Funding
We acknowledge the support by the Deutsche Forschungsgemeinschaft (DFG) under the grant Sh405/10 and Sh405/7 (SPP 1881 ‘Turbulent Superstructures’) and the Leibniz Supercomputing Centre (LRZ).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Theory for diffusion dominated flows
We follow the theory of Davey (Reference Davey1967), but solve the equations in a more general way, to allow for flexibility in the chosen BCs; for more details, the reader is referred to Davey (Reference Davey1967) or Kelly & Vreeman (Reference Kelly and Vreeman1970). Neglecting the mean vertical velocity component, assuming the mean horizontal velocity to be independent of 
$x$ and neglecting the contributions from the mean temperature field 
$\bar {\theta }$, the linearised, non-dimensionalised Navier–Stokes equations in two dimensions read
\begin{gather} \partial_t u^\prime + (U+u^\prime)\partial_x u^\prime + w^\prime \partial_z (U+u^\prime) ={-}\partial_x p + \nu^{*} \left( \frac{\partial^2 U}{\partial z^2} + \frac{\partial^2 u^\prime}{\partial x^2} + \frac{\partial^2 u^\prime}{\partial z^2}\right) , \end{gather}
\begin{gather}\partial_t w^\prime + (U+u^\prime)\partial_x w^\prime + w^\prime\partial_z (w^\prime) ={-}\partial_z p + \nu^* \left(\frac{\partial^2 w^\prime}{\partial x^2} + \frac{\partial^2 w^\prime}{\partial z^2}\right) + \theta^\prime , \end{gather}
\begin{gather}\partial_x u^\prime + \partial_z w^\prime = 0. \end{gather}
Here, 
$u^\prime$ and 
$w^\prime$ are, respectively, the horizontal and vertical components of the velocity fluctuations with respect to their time averages, i.e. 
$U$ and 
$W=0$, and 
$\theta ^\prime$ is the temperature fluctuation. For non-dimensionalisation we have used the free-fall velocity 
$u_{ff}\equiv (\alpha g \Delta \hat {H})^{1/2}$, the height 
$\hat {H}$ and the amplitude of the thermal TW, 
$\varDelta$, so that 
$\nu ^*=\sqrt {Pr/Ra}$. Let us consider a single wave mode in the horizontal 
$x$-direction and in time 
$t$, e.g.
\begin{gather} w^\prime(x,z,t) = \tfrac{1}{2}\left( \hat{w}(z)\exp({+\textrm{i}(kx-2{\rm \pi} \varOmega t)}) +\hat{w}^*(z) \exp({-\textrm{i}(kx-2{\rm \pi} \varOmega t)}) \right), \end{gather}
\begin{gather}u^\prime(x,z,t) ={-}\int \partial_z w^\prime\,{\textrm{d}x} = \frac{{\rm i}}{2k}\left( \partial_z\hat{w}(z) \exp({+\textrm{i}(kx-2{\rm \pi} \varOmega t)})\right.\nonumber\\ \quad-\left.\partial_z\hat{w}^*(z) \exp({-\textrm{i}(kx-2{\rm \pi} \varOmega t)}) \right), \end{gather}
\begin{gather}\theta^\prime(x,z,t) = \tfrac{1}{2}\left( \hat{\theta}(z)\exp({+\textrm{i}(kx-2{\rm \pi} \varOmega t)}) +\hat{\theta}^*(z)\exp({-\textrm{i}(kx-2{\rm \pi} \varOmega t)}) \right), \end{gather}where the asterisk denotes the complex conjugate of a function. We will consider two BCs (different scenarios), Scenario 1 describes a set-up where two travelling thermal waves are imposed at the top and the bottom (without any phase difference). This case was considered in the present work. Scenario 2, on the other hand, describes a set-up where the thermal wave travels only at the bottom, while the dimensionless top temperature equals zero.
Step 1: calculate 
$\hat {\theta }(z)$.
Neglecting dissipation in 
$x$, all convective terms and mean temperature contributions, the linearised non-dimensional energy equation reads
\begin{equation} \partial_t \theta^\prime = \kappa^*\left(\frac{\partial^2 \theta^\prime}{\partial z^2}\right), \end{equation}
where 
$\kappa ^*=1/\sqrt {RaPr}$. This, together with (A6), leads to the following equation for the wave amplitude equation 
$\hat {\theta }(z)$:
\begin{equation} \frac{\textrm{d}^2 \hat{\theta}}{\textrm{d}z^2} - \lambda^2 \hat{\theta} = 0; \quad \lambda^2 = \frac{2{\rm \pi} {\rm i}\varOmega}{\kappa^*}. \end{equation}The solution to (A8), for the two scenarios is
Scenario 1
\begin{align*} &\text{For}\ \hat{\theta}\vert_{z={-}1/2} = \hat{\theta}\vert_{z=1/2} = \tfrac{1}{2} :\\ &\hat{\theta}(z) = \frac{\cosh(\lambda z)}{2 \cosh(\lambda/2)} . \end{align*}
Scenario 2
\begin{align*} &\text{For}\ \hat{\theta}\vert_{z={-}1/2} = \tfrac{1}{2}, \hat{\theta}\vert_{z=1/2} = 0 :\\ &\hat{\theta}(z) = \frac{\sinh(\lambda/2 - \lambda z)}{2 \sinh(\lambda)} . \end{align*}
Step 2: calculate 
$\hat {w}(z)$.
Eliminate the pressure term by cross-differentiation of (A1) and (A2), substitute (A4)–(A6), neglect convective terms and assume that the thermal wavelength is much larger than the height of the cell (
$kH\ll 1$) to obtain
\begin{equation} \frac{\partial^4 \hat{w}}{\partial z^4} - \alpha^2 \frac{\partial^2 \hat{w}}{\partial z^2} = k^2\hat{\theta}, \quad \alpha^2 = \frac{2{\rm \pi} {\rm i}\varOmega}{\nu^*}. \end{equation}
For 
$\hat {w}\vert _{z=1/2} = \hat {w}\vert _{z=-1/2} = \partial _z \hat {w}\vert _{z=1/2} = \partial _z \hat {w}\vert _{z=-1/2} = 0$, the solution to (A9) is
\begin{equation} \hat{w}(z) = \frac{c_1}{\alpha^2}\cosh(\alpha z) +\frac{c_2}{\alpha^2}\sinh(\alpha z) + c_3 z + c_4 + c_5 \cosh(\lambda z) + c_6 \sinh(\lambda z). \end{equation}Scenario 1
\begin{align*} A &= \frac{k^2}{2\nu^* \lambda^2(\lambda^2-\alpha^2)}, \\ c_1 &={-} \lambda\alpha A \frac{\tanh(\lambda/2)}{\sinh(\alpha/2)}, \\ c_2 &= 0, \\ c_3 &= 0, \\ c_4 &= A\left( \frac{\lambda}{\alpha} \frac{\tanh(\lambda/2)}{\tanh(\alpha/2)} -1\right), \\ c_5 &= \frac{A}{\cosh(\lambda/2)}, \\ c_6 &= 0. \end{align*}
Scenario 2
\begin{align*} A &= \frac{k^2}{4\nu^* \lambda^2(\lambda^2-\alpha^2)}, \\ c_1 &={-} \lambda\alpha A \frac{\tanh(\lambda/2)}{\sinh(\alpha/2)}, \\ c_2 &= \frac{-\alpha A \left(\dfrac{\lambda}{\tanh(\lambda/2)} -2\right)}{(2/\alpha) \sinh(\alpha/2) - \cosh(\alpha/2)}, \\ c_3 &={-}\frac{c_2}{\alpha}\cosh(\alpha/2) + \frac{\lambda A}{\tanh(\lambda/2)}, \\ c_4 &= A\left( \frac{\lambda}{\alpha} \frac{\tanh(\lambda/2)}{\tanh(\alpha/2)} -1\right), \\ c_5 &= \frac{A}{\cosh(\lambda/2)}, \\ c_6 &= \frac{-A}{\sinh(\lambda/2)}. \end{align*}
Step 3: calculate 
$U(z)$.
Averaging equation (A1) over time and over one wavelength in 
$x$, we obtain the following equation for the mean flow 
$U(z)$:
\begin{equation} \nu^* \frac{\textrm{d}^2 U}{\textrm{d}z^2} = \frac{\textrm{d}}{\textrm{d}z} (\overline{u^\prime w^\prime}), \end{equation}which can be solved via numerical integration using the no-slip BCs at the plates.
In addition, in the supplementary material we provide a Python code snippet, which gives the solution for the various quantities 
$\hat {\theta }, \hat {u},\hat {w},\langle u^\prime w^\prime \rangle$. Note that 
$z$ runs from 
$-1/2$ to 
$1/2$ and there is a singularity for 
$Pr= 1$, which can be avoided by choosing a value very close to one or could be resolved by L'Hôpital's rule.
Appendix B. Heat and momentum transport
The Nusselt number 
$Nu$ and Reynolds number 
$Re$, based on the wind velocity, are defined as
\begin{equation} Nu \equiv{-} \left\langle\frac{\partial \bar{\theta}}{\partial z}\biggr\rvert_{z=0} \right\rangle_{\mathcal{A}}, \quad Re \equiv \sqrt{\frac{Ra}{ Pr}} \sqrt{\langle \overline{\boldsymbol{u}^2}\rangle_{V}}, \end{equation}
where 
$\mathcal {A}$ denotes the horizontal plane for the cylinder or, respectively, the 
$x$-direction for the 2-D simulations. Figure 13 shows 
$Nu(\varOmega )$ and 
$Re(\varOmega )$, normalised by their values at 
$\varOmega =10^{-3}$. Their exact values are given in the supplementary material. The 2-D system (figure 13a,b) shows a significant heat and momentum transport enhancement for certain TW speeds 
$\varOmega$, especially for large 
$Ra$. For the 3-D cylindrical system (figure 13c), no clear correlation between the zonal flow maximum (see figure 9) and 
$Nu(\varOmega )$ and 
$Re(\varOmega )$ is observed. However, a small 
$Re$ enhancement is present at 
$\varOmega \approx 10^{-2}$.
Figure 13. Normalised 
$Nu$ and 
$Re$ vs. 
$\varOmega$ for (a) 2-D set-up A, (b) 2-D set-up B and (c) 3-D cylinder; 
$Ra =10^3$ (
$\bullet$, blue), 
$10^4$ (
$\bullet$, orange), 
$10^5$ (
$\bullet$, green), 
$10^6$ (
$\bullet$, red) and 
$10^7$ (
$\bullet$, black).
Appendix C. Linear stability analysis
In § 3.1.1 a temporal linear stability analysis was conducted to identify the most unstable eigenmode of the 2-D linearised Navier–Stokes equations, where a wave-like form was considered only in time. Thus, any flow quantity 
$\phi (x,z,t)$ is represented as 
$\phi (x,z,t) = \hat {\phi }(x,z) \textrm {e}^{-\textrm {i}\omega t}$ and the system of equations for the horizontal velocity 
$u$, the vertical velocity 
$w$, the pressure 
$p$ and the temperature 
$\theta$ reads
\begin{equation} \begin{bmatrix} L_{2D} + D_x U & D_z U & D_x & 0\\ D_x W & L_{2D} + D_z W & D_z & -1\\ D_x & D_z & 0 & 0\\ D_x \bar{\theta} & D_z \bar{\theta} & 0 & K_{2D}\\ \end{bmatrix} \begin{bmatrix} \hat{u} \\ \hat{v} \\ \hat{p} \\ \hat{\theta} \end{bmatrix} = \omega \begin{bmatrix} \textrm{i} & 0 & 0 & 0 \\ 0 & \textrm{i} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \textrm{i} \\ \end{bmatrix} \begin{bmatrix} \hat{u} \\ \hat{v} \\ \hat{p} \\ \hat{\theta} \end{bmatrix}, \end{equation}where
\begin{gather} L_{2D} = U D_x + W D_z + \sqrt{Pr/Ra}\left({-}D_x^2-D_z^2\right) , \end{gather}
\begin{gather}K_{2D} = U D_x + W D_z + 1/\sqrt{RaPr}\left({-}D_x^2-D_z^2\right). \end{gather}The overline represents the mean field quantity. In our study we applied the Chebyshev method to approximate the vertical gradient (
$D_z$) and the Fourier method for the horizontal gradient (
$D_x$). Conveniently, the corresponding differentiation matrices are available open source, e.g. we used the Python package dmsuite.
The linear set of (C1) is solved as a generalised eigenvalue problem of the form 
$\boldsymbol {A} \hat {\phi }=\omega \boldsymbol {B} \hat {\phi }$, where the eigenvectors 
$\phi (x,z,t)$ represent the wave amplitudes and the eigenvalues 
$\omega$ their respective temporal behaviour. The matrices, of size 
$4 \times N_x \times N_z$, are very large and therefore an iterative solver has to be used (e.g. Python's scipy.eigs). The code has been validated by solving the Blasius boundary layer, pipe flow and Rayleigh–Taylor instabilities in one and two dimensions, and in closed and periodic domains. For all cases we have found excellent agreement with results in the literature.
Open access
