From zonal flow to convection rolls in Rayleigh-B\'enard convection with free-slip plates

Rayleigh-B\'enard (RB) convection with free-slip plates and horizontally periodic boundary conditions is investigated using direct numerical simulations. Two configurations are considered, one is two-dimension (2D) RB convection and the other one three-dimension (3D) RB convection with a rotating axis parallel to the plate. We explore the parameter range of Rayleigh numbers Ra from $10^7 to $10^9$ and Prandtl numbers $Pr$ from $1$ to $100$. We show that zonal flow, which was observed, for example, by Goluskin \emph{et al}. \emph{J. Fluid. Mech.} 759, 360-385 (2014) for $\Gamma=2$, is only stable when $\Gamma$ is smaller than a critical value, which depends on $Ra$ and $Pr$. With increasing $\Gamma$, we find a second regime in which both zonal flow and different convection roll states can be statistically stable. For even larger $\Gamma$, in a third regime, only convection roll states are statistically stable and zonal flow is not sustained. For the 3D simulations, we fix $Ra=10^7$ and $Pr=0.71$, and compare the flow for $\Gamma=8$ and $\Gamma = 16$. We demonstrate that with increasing aspect ratio $\Gamma$, zonal flow, which was observed for small $\Gamma=2\pi$ by von Hardenberg \emph{et al}. \emph{Phys. Rev. Lett.} 15, 134501 (2015), completely disappears for $\Gamma=16$. For such large $\Gamma$ only convection roll states are statistically stable. In between, here for medium aspect ratio $\Gamma = 8$, the convection roll state and the zonal flow state are both statistically stable. What state is taken depends on the initial conditions, similarly as we found for the 2D case.

Rayleigh-Bénard (RB) convection with free-slip plates and horizontally periodic boundary conditions is investigated using direct numerical simulations. Two configurations are considered, one is two-dimension (2D) RB convection and the other one three-dimension (3D) RB convection with a rotating axis parallel to the plate, which for strong rotation mimics 2D RB convection.
For the 2D simulations, we explore the parameter range of Rayleigh numbers Ra from 10 7 to 10 9 and Prandtl numbers P r from 1 to 100. The effect of the width-to-height aspect ratio Γ is investigated for 1 Γ 128. We show that zonal flow, which was observed, for example, by Goluskin et al. J. Fluid. Mech. 759, 360-385 (2014) for Γ = 2, is only stable when Γ is smaller than a critical value, which depends on Ra and P r. With increasing Γ , we find a second regime in which both zonal flow and different convection roll states can be statistically stable. For even larger Γ , in a third regime, only convection roll states are statistically stable and zonal flow is not sustained. How many convection rolls form (or in other words, what the mean aspect ratio of an individual roll is), depends on the initial conditions and on Ra and P r. For instance, for Ra = 10 8 and P r = 10, the aspect ratio Γ r of an individual, statistically stable convection roll can vary in a large range between 16/11 and 64. A convection roll with an as large aspect ratio of Γ r = 64 or more generally already with Γ r 10 can be seen as localized zonal flow, and indeed carries over various properties of the global zonal flow.
For the 3D simulations, we fix Ra = 10 7 and P r = 0.71, and compare the flow for Γ = 8 and Γ = 16. We first show that with increasing rotation rate both the flow structures and global quantities like the Nusselt number N u and the Reynolds number Re increasingly behave like in the 2D case. We then demonstrate that with increasing aspect ratio Γ , zonal flow, which was observed for small Γ = 2π by von Hardenberg et al. Phys. Rev. Lett. 15, 134501 (2015), completely disappears for Γ = 16. For such large Γ only convection roll states are statistically stable. In between, here for medium aspect ratio Γ = 8, the convection roll state and the zonal flow state are both statistically stable.

Introduction
Large scale so-called zonal flows, which display strong horizontal winds, can be observed in many buoyancy-driven flows. Typical examples include zonal flow in the atmosphere of Jupiter (Heimpel et al. 2005;Kaspi et al. 2018;Kong et al. 2018) and other three Jovian planets (Ingersoll 1990;Sun et al. 1993;Cho & Polvani 1996;Yano et al. 2003), in the oceans (Maximenko et al. 2005;Richards et al. 2006;Nadiga 2006), and possibly in the Earth's outer core (Miyagoshi et al. 2010). In toroidal tokamak devices, zonal flows in the poloidal direction are crucial to confine plasmas magnetically (Diamond et al. 2005).
How to study such flows? In general, Rayleigh-Bénard (RB) convection Chillà & Schumacher 2012), i.e., a fluid in a container heated from below and cooled from above, is the paradigmatic model system for buoyancydriven flows. The key control parameters are the Rayleigh number Ra = gα∆H 3 /νκ and the Prandtl number P r = ν/κ. Here, g is the gravitational acceleration, α the thermal expansion coefficient, H the height of the fluid sample, ∆ = T b − T t the temperature difference between the hot bottom and the cold top plate, κ the thermal diffusivity, and ν the kinematic viscosity of the fluid. The third control parameter is the aspect ratio Γ , which is defined as the ratio of the width to the height of the container. The response of the system is characterized by the Nusselt number N u = QH/(k∆) and the Reynolds number Re = U H/ν, which indicate the non-dimensional heat transport and flow strength in the system, respectively. Here Q is the heat flux crossing the system and U = u · u V,t the characteristic velocity, where V,t indicates volume and time averaging. Indeed, to study zonal flow, RB convection with free-slip plates and horizontally periodic boundary conditions has commonly served as model system (Goluskin et al. 2014;van der Poel et al. 2014;von Hardenberg et al. 2015;Novi et al. 2019). Here the free-slip at the plates is very crucial to enable the zonal flow; for no-slip boundary conditions, zonal flow is significantly suppressed and it only exists for very small Γ (van der Poel et al. 2014).
In the 2D version of the RB system with free-slip plates and horizontally periodic boundary conditions, indeed, for small Γ = 2, zonal flow develops readily since the horizontal periodicity allows for a horizontal mean flow, while the free-slip boundaries apply no shear stresses to slow down the fluid. In addition, the two-dimensionality precludes transverse perturbations that could disrupt the mean flow (van der Poel et al. 2014;Goluskin et al. 2014). Such zonal flow in 2D RB convection has attracted quite some attention because of its relevance to thermal convection in the atmosphere (Seychelles et al. 2008(Seychelles et al. , 2010Bouchet & Venaille 2012). For free-slip boundary conditions at the plates, P r = 1, Ra 10 7 , and a small Γ = 2, van der Poel et al. (2014) found that a flow topology consisting of two shear layers with a predominant horizontal motion is formed. The flow in the lower half of the domain moves in the opposite direction as that in the top part. Most of the time, the heat transfer of this flow is N u ≈ 1, while there are intermittent bursts in which N u 1. Goluskin et al. (2014) studied the characteristics of such 2D zonal flows in a periodic Γ = 2 cell for an extended parameter range 10 3 Ra 10 10 and 1 P r 10. They found that for P r 2, the zonal flow undergoes strong global oscillations on long timescales. Also intermittent bursts in the heat transport as in van der Poel et al. (2014) are observed. For P r 3, the zonal flow is sustained at all times without bursts, and the Nusselt number N u is always much larger than 1.
To what degree can 2D simulations mimic the dynamics in three-dimensional (3D) flows? Actually many 3D geophysical and astrophysical flows exhibit certain 2D properties when anisotropic effects, such as geometrical confinement, rapid rotation, stratification, or magnetic fields, are imposed. We will show in this study how the 2D flow structures arise with increasing rotation rate for RB convection rotating about a horizontal axis. Such flow will be called span-wise rotating RB convection in this paper. 2D simulations, which are computationally more accessible than 3D simulations, have also been widely used to study RB convection with no-slip plates (Johnston & Doering 2009;Sugiyama et al. 2010;Huang & Zhou 2013;van der Poel et al. 2015b;Zhu et al. 2018aZhu et al. , 2019Wan et al. 2020). In van der Poel et al. (2013) 2D and 3D simulations are compared in detail, and many similarities are found for P r 1.
In contrast to in 2D configuration, zonal flow has not been reported in horizontally isotropic 3D simulations of RB convection with free-slip plates (Petschel et al. 2013;Kunnen et al. 2016). It seems that in 3D convection zonal flow only appears when an anisotropy is added. For example, von Hardenberg et al. (2015) found that a strong zero-wavenumber wind (i.e., zonal flow) can arise in 3D RB convection if the horizontal isotropy is broken by strong enough uniform rotation about a horizontal axis. Indeed, according to the Taylor-Proudman theorem, the flow can become 2D-like when the rotation is sufficiently fast. Recently, Novi et al. (2019) further generalized the situation by varying the tilting angle of the rotation axis with respect to gravity. This configuration mimics the flow at different latitudes in a rotating fluid shell. A large-scale cyclonic vortex tilted along the rotation axis is identified for φ between 45 • and 90 • , where φ is the angle between the rotation axis and the horizontal plane. At moderate latitudes the calculations of Novi et al. (2019) suggest the possible coexistence of zonal jets and tilted-vortex solutions.
Even though flows in geophysics, astrophysics, and plasma-physics often occur in largeaspect ratio systems, most of the previous zonal flow studies with free-slip conditions at the plates were performed for horizontally periodic small-aspect ratio cells, typically Γ = 2 (2D) or Γ = 2π (3D). However, recent studies on large-aspect ratio 3D RB convection with no-slip plates revealed the existence of superstructures with horizontal extent of six to seven times the height of the domain (Hartlep et al. 2003;Parodi et al. 2004;Pandey et al. 2018;Stevens et al. 2018;Krug et al. 2020;Green et al. 2020). These findings motivated us to study zonal flow at much larger Γ (up to 128) than hitherto done, in order to test whether zonal flow will sustain at these much larger Γ , or whether some other large-scale structures evolve, which are not captured in simulations with small Γ .
We will find that for free-slip plates and periodic boundary conditions the aspect ratio indeed has a very strong influence on the flow phenomena in 2D RB convection and in 3D RB convection with span-wise rotation. In particular, we will show that zonal flow is only stable when the aspect ratio of the system is smaller than a critical value, which depends on Ra and P r; it disappears in large-aspect ratio flow configurations.
The paper is organized as follows. We first describe the simulation details in Section 2. The 2D results are presented and analysed in Section 3, which is divided into three parts. Section 3.1 demonstrates the disappearance of zonal flow with increasing the aspect ratio Γ . Section 3.2 studies the coexistence of multiple convection roll states. The effective scaling relations for N u(Ra, P r) and Re(Ra, P r) for different convection roll states are discussed in Section 3.3. The 3D RB convection with increasing rotation strength about an axis parallel to the plate (i.e., increasing two-dimensionalization according to Taylor Table 1: Overview of the run simulations in 2D. The first three columns indicate the Ra, P r, and Γ range of the simulations. N r z and N z z indicate the number of grid points in the vertical direction for the simulations with initial conditions of roll states and shear flow, respectively. N r BL and N z BL indicate the minimum number of grid points in the thermal boundary layer for convection roll states and zonal flow state, respectively. We note that the number of grid points in the boundary layer is always higher than that given by the recommendation of Shishkina et al. (2010) for the no-slip case, which is about 5 to 9 for this Ra range, to ensure that the boundary layers are resolved. The number of grid points in the horizontal direction is generally equal or larger than N x = N z × Γ . For Ra = 10 8 and 3 × 10 8 , N z z = 256 is used only for the large Γ cases (fore example, Ra = 3 × 10 8 , Γ = 32 and 64) where very long simulations are performed, in order to test whether the zonal flow state can stably exist. Figure 1: Sketch of (a) 2D RB convection and (b) 3D RB convection with span-wise rotation for free-slip plates and horizontally periodical conditions. Proudman theorem) is discussed in Section 4, where we also show the transition from zonal flow to convection roll states with increasing Γ . We summarize our findings in Section 5.

Simulation details
The configurations and the coordinate systems used in this work are shown in figure 1. We performed direct numerical simulations using the second-order staggered finite difference code AFiD. Details about the numerical method can be found in (Verzicco & Orlandi 1996;van der Poel et al. 2015a;Zhu et al. 2018b). The governing equations in dimensionless form read From zonal flow to convection rolls in RB convection with free-slip plates 5 ∇ · u = 0, (2.1) Here e y and e z are the unit vector in the y and z direction, respectively. u, t, p, θ are velocity, time, pressure and temperature, respectively. The length and velocity are nondimensionalized using the height of the convection cell H and the free-fall velocity U = (gα∆H) 1/2 , respectively. This implies as reference time the free-fall time t f = H/U . For the 3D simulations also the Rossby number Ro = U/(2ΩH) is used, where Ω is the angular velocity. Non-uniform grids with points clustered near the top and bottom plates are employed.
How to choose the initial conditions to trigger the different flow states? For the zonal flow simulations we used a linear shear-flow velocity profile u(z) = 2z − 1, w = 0 in combination with a linear temperature profile θ(z) = 1 − z with additional random perturbations as initial conditions. In addition, different convection roll states were generated using a Fourier basis: u(x, z) = sin(n (i) πx/Γ )cos(πz), w(x, z) = −cos(n (i) πx/Γ )sin(πz) where n (i) indicates the number of initial rolls in the horizontal direction, while the initial temperature is the same as zonal flow simulations. A similar Fourier basis was also used to study heat transport (Chong et al. 2018) and flow reversals (Chandra & Verma 2011;Wang et al. 2018bWang et al. , 2019aChen et al. 2019). An overview of the 2D simulations and the used grid resolutions are given in table 1. Table 3 and 4 in the Appendix list 2D simulation details for the main cases where N u and Re are discussed. The 3D simulation details are given in table 5, where the 2D simulations for the corresponding parameters are also listed for comparison.

Disappearance of zonal flow with increasing Γ
We first show what will happen to zonal flow with increasing Γ . From Goluskin et al. (2014), it is known that zonal flow exists for Ra = 10 8 , P r = 10, Γ = 2. With increasing Γ , we find that for this Ra and P r, zonal flow can stably exist for Γ 12. In figure 2(a) we show that zonal flow is statistically stable at least up to 200 000 free-fall time units for Γ = 4 and at least up to 100 000 free-fall time units for Γ = 12. Here, we used "statistically stable" to denote that the corresponding chaotic flow state is always sustained in our long-time simulations. Temperature snapshot of the zonal flow for Γ = 12 in figure 2(b) demonstrates that hot plumes drift leftwards, and cold plumes drift rightwards. This produces a strong horizontal shear in which however the vertical heat transport is low.
We now explore even larger aspect ratio domains. Figure 2(c) shows the time evolution of N u for three separate simulations with random perturbations added to the initial temperature field for Ra = 10 8 and P r = 10 in a Γ = 64 cell. For all the three simulations, the zonal flow eventually evolves to a convection roll state. The time at which the transition occurs is very different for each simulation. The reason is that the flow is susceptible to small differences in the initial conditions, which are different for each simulation due to the random perturbations to the initial temperature field. Such  The three curves correspond to three separate simulations with random perturbations added to the initial temperature field. In all the cases, the flow undergoes a transition from zonal flow to convection roll states, for which N u is larger. (d ) Temperature snapshots at different times denoted by the red dashed lines for the simulation indicated by the red curve in panel (c). At t = 2000, there is zonal flow, whereas later it features an increasing number of turbulent convection rolls. (e) The final two-roll state for Ra = 10 8 , P r = 10, Γ = 128, and the zoom in of the two plumeejecting regions. For all these simulations the initial velocity had a linear shear flow profile u(z) = 2z − 1, w = 0, in order to trigger a zonal flow state. sensitivity to the initial conditions is typical for chaotic systems and makes it impossible to predict when the transition happens. Figure 2(d) displays three temperature snapshots at different time instants indicated by the red dashed lines in figure 2(c). A complementary movie, showing how the zonal flow undergoes a transition towards a convection roll state, is given in the supplementary material. Initially, the hot plumes travel leftwards and the cold plumes to the right. The transition starts when some local hot plumes are strong enough to deviate upwards and cross the whole fluid layer up to the collision with the upper cold plate. This prevents the further rightward motion of the neighbouring cold plumes, which instead start to move downwards. This process generates a local large-scale circulation, as is observed in the temperature field at t = 4172. The circulation grows over time until two stable convection rolls of equal size are formed as seen in the temperature filed at t = 5000. Figure 2(e) shows that we also obtain a two-roll state for Γ = 128, such that the horizontal extent of the convection roll is 64 times the height of the convection cell, and this state is stable for more than 10 000 free-fall time units.
We now explore the phase diagram for the different flow states in the parameter space spanned by Ra, P r, and Γ , see figure 3 for the simulated cases. We find that in smallaspect ratio cells, only zonal flow is stable, while in large-aspect ratio cells, only convection roll states are stable. For intermediate aspect ratios, we find a regime in which both zonal flow and convection roll states are stable, depending on the initial conditions mentioned in Section 2. We call this regime the bi-stable one. In order to map out an accurate phase diagram, we performed long simulations for the bi-stable cases with largest Γ to conclude that the corresponding zonal flow can stably exist and does not evolve to a convection roll state. Overall, we ran the simulations at least 50 000 free-fall time units for these cases, which corresponds to at least 5 viscous diffusive time units (H 2 /ν) or 0.5 thermal diffusive time units (H 2 /κ). From figure 3(a) it can be seen that, when Ra is increased, the bi-stable regime exists in an increasing Γ range. This is also consistent with the finding that zonal flow develops more readily for higher Ra for Γ = 2 (Goluskin et al. 2014). Figure 3(b) demonstrates that the Γ range for the bi-stable state also depends on P r. For Ra = 10 8 the largest range of bi-stable state exists at P r ∼ 30, namely between Γ = 3 and Γ = 24.
We have already shown that zonal flow can not be sustained, and only convection roll states are observed, when Γ is larger than a critical value, which depends on Ra and P r. A related question is how many convection rolls (in other words, what is the mean aspect ratio of individual convection roll) can develop for a specific (Ra, P r, Γ ). In the next subsection, we will explore the possible convection roll states using different initial roll states generated by different Fourier basis, as explained in Section 2.

Coexistence of multiple convection roll states at large Γ
In this subsection, we study the coexistence of multiple convection roll states in largeaspect ratio domains, all being statistically stable states once achieved. Figure 4 shows that for Ra = 10 8 , P r = 10, and free-slip at the plates, in a Γ = 16 system convection rolls with a mean dimensionless horizontal size of 1.6 Γ r 8 are all statistically stable. The heat transport considerably increases with decreasing mean aspect ratio Γ r of an individual convection roll. For example, the heat transport for the Γ r = 1.6 state is almost twice as high as that for Γ r = 8. Although it had been observed before that 8 Q. Wang et al. Figure 4: Temperature snapshots of different roll states for Ra = 10 8 and P r = 10 in a Γ = 16 periodic cell. The dimensionless mean horizontal size of the convection roll Γ r (i.e., the mean aspect ratio of one individual roll) and the Nusselt number N u for each state are indicated. The different roll states are from initial conditions with different number of initial rolls. All these states can stably exist for long time (see table 3 in Appendix) without undergoing a transition to other states.
convection rolls with smaller Γ r imply a higher heat transport -e.g. for 2D RB convection with no-slip plates (van der Poel et al. 2012;Wang et al. 2018a), for RB convection in an annulus convection cell (Xie et al. 2018), and also for Taylor-Couette flow (Huisman et al. 2014) -in those cases the observed increase in the transport is typically a couple of percents, and by far not as large as observed for RB with free-slip plates and large aspect ratio cells as studied here. This difference is due to different plume dynamics and the associated spatial dependence of the local Nusselt number, N u(x) , as we will discuss later.
We also tested initial conditions with 12, 14, 16 rolls for Ra = 10 8 , P r = 10, Γ = 16. These states with smaller rolls are not stable and will finally undergo a transition to the ten-roll state by merging of convection rolls. From figure 5(a) is is seen that the vertical Reynolds number Re z has a sudden decrease during merging of rolls, because the strong vertical motion is concentrated near the plume-ejecting regions between two neighbouring rolls. The decrease of N u during merging events can also be observed in figure 5(b), which is related to the decreased vertical motion. Figure 5(c) shows how the flow undergoes a transition from the initial sixteen-roll state to the final ten-roll state by successive merging of convection rolls. The transition happens when the balance of the roll state is broken by horizontal motion of local hot/cold plumes: In the second snapshot at t = 85 the system is still in the initial sixteen role state. However, one can already see that a local hot plume moves leftwards while its neighbouring cold plume moves rightwards (marked by red arrows). In the third snapshot at t = 87 two hot plumes merge to a single one and so do two cold ones, thus annihilating two rolls. The resulting fourteen-roll state is shown in the fourth snapshot taken at t = 141. At later times the horizontal motion of the plumes and their further merging let the fourteen-roll state evolve to a twelve-roll state (t = 182, III), and finally to a ten-roll state (t = 300, IV). We also performed very long simulations as indicated in figure 5(d), from which we conclude that the ten-roll state can statistically stably exist for a very long time without undergoing any further transition to yet another state. Figure 6 displays phase diagrams for all the possible convection roll states in the Ra−Γ r is horizontal Reynolds number and Re z = (Ra/P r) ( w 2 V ) the vertical one. (c) Temperature snapshots at different times. The roll merging can be seen, namely the flow undergoes a transition from initial sixteen-roll state (I), to a fourteen-roll state (II), to a twelve-roll state (III), and then to the final ten-roll state (IV). The figure has the same colour scale as figure 4. (d ) Time evolution of Re for much longer time (on a log-scale) to show that the final ten-roll state is stable without undergoing a transition to an other roll state. Circles denote that the corresponding roll state with the mean aspect ratio Γ r of an individual roll is stable, while crosses denote that the roll state is not stable. The solid line in (b) connects the minimal mean aspect ratio Γ r,min of an individual convection roll for different P r for Γ = 16, while the dashed line connects Γ r,min for different P r for Γ = 32. and P r−Γ r parameter spaces. The stable roll states can last for several thousand free-fall time units without undergoing a transition to other states (see Appendix). Figure 6(a) shows a weak dependence of Γ r on Ra. One can observe the same stable roll states for the considered Ra range. In contrast, a pronounced dependence of Γ r on P r is observed in figure 6(b), where convection roll states with the smallest Γ r are observed for intermediate P r ≈ 10. The minimum Γ r = 16/11 occurs for Γ = 32, which means that the horizontal extent of a stable convection roll is always larger than the height of the system. This explains why convection rolls cannot be supported for small Γ ≈ 2, where indeed only zonal flow was obtained. From figure 6(b) it can also be concluded that these results are independent of the aspect ratio of the system once it is large enough, as we obtained almost the same result for Γ = 16 and Γ = 32 domains.

Nusselt number and Reynolds number
We now discuss the effective scaling relations of N u and Re as function of Ra and P r. These relationships are usually expressed with effective scaling laws N u ∼ Ra γ N u P r α N u and Re ∼ Ra γ Re P r α Re ). The effective scaling laws have been widely discussed for no-slip cases for both 2D and 3D convection ). For the 2D horizontally periodic cases with no-slip plates, N u ∼ Ra 0.29 is found with P r = 1, Ra 10 10 (Johnston & Doering 2009;Zhu et al. 2018a  and Re ∼ Ra 0.6 (Sugiyama et al. 2009;Zhang et al. 2017;Wang et al. 2019b). However, how these effective scaling relations will change for free-slip plates has not been explored, especially not for convection roll states, which are only present in large enough domain size. Figure 7(a) and 7(b) show N u and Re as function of Ra for both zonal flow (Γ = 2) and convection roll states (Γ = 16) for P r = 10. Figure 7(a) reveals that the heat transfer in the convection roll state is much higher than that of the zonal flow state. Detailed information about the obtained scaling exponents is listed in table 2. For the convection roll states we find that the effective scaling exponent γ N u in N u ∼ Ra γ N u is approximately 0.3. It increases with increasing Γ r , reaching about 1/3 for the largest Γ r = 8, which is the value predicated by the Grossmann-Lohse (GL) theory for the I < ∞ and III ∞ regimes (Grossmann & Lohse 2000Shishkina et al. 2017) for the no-slip case. For zonal flow γ N u is much smaller, namely only 0.17. The effective scaling exponent γ Re in Re ∼ Ra γ Re is about 0.6 for zonal flow and about 0.67 for the convection roll state, which is also close to the GL predication of 2/3 for the I < ∞ and III ∞ regimes (Grossmann & Lohse 2001;Shishkina et al. 2017) for no-slip cases, while it is larger than 0.6 that has been reported for 2D RB convection with no-slip plates (Sugiyama et al. 2009;Zhang et al. 2017;Wang et al. 2019b).
Next, we will discuss the Prandtl number dependence of the Nusselt number, N u(P r). For no-slip plates in 3D cases, it is known that the N u is maximal around P r ∼ 2 − 3, and after that it declines with increasing P r (Ahlers & Xu 2001;Xia et al. 2002;Stevens et al. 2011). This remarkable maximum in N u(P r) had in fact been predicted before by the GL theory (Grossmann & Lohse 2000. In contrast, for the 2D cases, Huang & Zhou (2013) showed that N u(P r) has a minimum, rather than maximum as in the 3D case, namely at P r ∼ 2 − 3 for moderate Ra. This anomalous relation is caused by counter-gradient heat transport in 2D cases.
How does the N u(P r) dependence look like for the 2D RB case with free-slip plates? For the zonal flow state Goluskin et al. (2014) already showed that N u is an increasing function of P r in the range 1 P r 10. Figure 7(c) and 7(d) show the relations for N u(P r) and Re(P r) for all states with free-slip plates (i.e., both for zonal flow and for various convection roll states), where we included both Γ = 16 and 32 data for the convection roll state. From figure 7(c) it can be seen that the N u(P r) trend shown by Goluskin et al. (2014) is also valid for the wider range of P r analysed in this paper. The reason why N u is much smaller for small P r is that zonal flow features intermittent 12 Q. Wang et al. bursts whereas most of the time N u is around 1 (Goluskin et al. 2014). For large P r, the flow does not burst and convective heat transport with N u 1 is sustained at all times, thus the corresponding N u is larger than that of the small P r cases.
For the convection roll states, figures 7(c) and 7(d) reveal that both N u and Re collapse well for the same mean aspect ratio of an individual convection cell (Γ r = 8, 4, 8/3) for both Γ = 16 and 32. This suggests that the global transport properties are directly related to the mean convection roll size Γ r . Figure 7(c) also shows that N u increases monotonically with increasing P r for large mean convection roll size Γ r = 16. This can be interpreted as that the flow for the large Γ r = 16 cases can be viewed as "localised" zonal flow, thus the N u(P r) follows the similar trend as that of zonal flow. In contrast, for small Γ r , N u decreases with increasing P r. The different N u(P r) trends for large and small Γ r can also be seen in the inset of figure 7(c). For Re(P r), figure 7(d) shows that the Re follows Re ∼ P r −1 for convection roll states (see table 2), the exponent −1 is the same as that of the GL predication for the I < ∞ and III ∞ regimes for no-slip cases. In order to understand different N u(P r) dependence for large and small Γ r as shown in figure 7(c), we first look at the flow organizations for convection roll states for different P r. Figure 8(a) gives the time-averaged temperature fields for the Γ r = 16 roll state for different P r. The flow near the bottom plate can be divided into plume-ejecting region, plume-impacting region, and between them there is wind-shearing region which occupies large fraction of the domain. In the ejecting region, thermal plumes are emitted, while in the impacting region, the boundary layer is impinged by the plumes from the opposite plate. The wind-shearing region is sheared by the large-scale circulation. The impacting regions on the top plate are the opposite of the ejecting regions of the bottom plate and vice versa. This kind of division is also adopted in periodic 2D RB convection with no-slip plates (van der Poel et al. 2015b). A remarkable observation is the stable-stratification near the plume impacting region. Figure 8(c) shows a zoom-in of the regions where hot plumes are ejected for P r = 1 and 10. When hot fluid impinges the cold plate, it does not have sufficient time to cool down, before it moves horizontally. The consequence is that the temperature of the fluid between the top boundary layer and the bulk is higher than that of the bulk fluid, thus implying a stable-stratification. This behavior is even observed at the centre line of the hot plume (the vertical line at the horizontal location where the local bottom wall heat flux is minimal) as is shown in temperature profiles in figure 8(d), where for P r > 1 stable-stratification near the cold plate is observed. The stable-stratification has also been observed at the axis in cylindrical RB convection (Tilgner et al. 1993;Brown & Ahlers 2007;Wan et al. 2019), and in 2D RB convection with no-slip plates and sidewalls for unit aspect ratio in the central region near the plates (Wan et al. 2020). The instantaneous temperature fields shown in figure 8(b) for P r = 100 reveals the "localised" zonal flow structures. It can be seen that plumes are ejected everywhere while they can only move vertically and impinge the cold plate in the central region.
Next, we focus on the local properties of the wall heat flux. The local wall heat flux is expressed through the local wall Nusselt number N u(x)| z=0,1 = ∂ θ t /∂z| z=0,1 . Figure  8(e) and 8(f ) show the spatial N u(x) dependence at the plates for Ra = 10 8 and Γ = 32 for Γ r = 16 and Γ r = 4, respectively. For small P r = 1, one sees from figure 8(a) the accumulation of hot fluid in the ejecting region near the bottom plate, which cause small temperature gradient (see figure 8(d)), and correspondingly small local N u (see figure  8(e)). So the centre of the ejecting region can be denoted as the point where local wall heat flux is minimal. In contrast, for the impacting region (x/H ≈ 0 at the bottom plate) where cold fluid directly impinges the hot plate, there is a sharp temperature gradient and thus large local N u (see figure 8(e)). Similar N u behavior is also observed in the Panels (e) and (f ) show the spatial dependence of N u(x) at the hot plate at z = 0 (solid lines) and the cold plate at z = 1 (dashed lines) for different P r for the (e) Γ r = 16 and the (f ) Γ r = 4 roll states. Note that all curves are shifted such that the minimum local N u at the hot plate is located at x/H = 16. ejecting/impacting regions near the top plate (dashed lines in figure 8(e)). The physical interpretation is as follows: The heat is ejected into the system through the bottom plate mainly at the plume-impacting regions where local temperature gradient is large, and then it is advected by large-scale circulation to the plume-ejecting regions, where conductive heat flux is low on the wall while convective heat flux is high above the wall. The heat is mainly removed from the system when hot plume impinges the cold plate.
For Γ r = 16, there is only one impacting region near the bottom plate, and the heat input is still dominated by the wind-shearing region, which occupies a large fraction of the domain. As the wind-shearing region is like "localised" zonal flow where N u increases with increasing P r, the global N u thus also increases with increasing P r. In contrast, for smaller Γ r , there are more impacting regions on the bottom plate, and these impacting regions contribute significantly to the global heat input. As the heat flux at the impacting region increases with decreasing P r, we thus see that the global N u also increases with decreasing P r.

Three-dimensional simulations
We have already shown that the zonal flow observed in 2D RB convection with freeslip plates and horizontally periodic boundary conditions for Γ = 2 (van der Poel et al. 2014;Goluskin et al. 2014) eventually disappears with increasing Γ . What about in 3D, under the same conditions? For the 3D RB convection with free-slip plates, previous studies have not reported zonal flow (Petschel et al. 2013;Kunnen et al. 2016). However, if we introduce span-wise rotation as illustrated in figure 1(b) where the rotating axis is parallel to the y axis, the flow will become 2D-like at sufficiently fast rotation, due to the Taylor-Proudman theorem. In this way we may observe zonal flow at certain parameters, as indeed was already reported in (von Hardenberg et al. 2015). Therefore here we will study span-wise rotating RB convection, focusing on the transition from zonal flow to convection roll states with increasing aspect ratio Γ of the container.
We first show that both the global transport properties like the Nusselt number N u and the Reynolds number Re, as well as the flow organization increasingly behave like 2D cases with increasing rotation rate. We fix the Rayleigh number to Ra = 10 7 and the Prandtl number to P r = 0.71. To be on the safe side, we choose a large domain with Γ = 16, as previous studies showed that large aspect ratios are needed in order to capture the superstructures which have a horizontal size of 6-7 times the height of the domain for 3D RB convection with no-slip plates (Pandey et al. 2018;Stevens et al. 2018;Krug et al. 2020;Green et al. 2020). The initial conditions have zero velocity and a linear temperature profile for these simulations. For RB convection rotating about a vertical axis for small P r = 0.71 with no-slip plates, N u initially does not change much with increasing the rotation rate (denoted by the inverse Rossby number 1/Ro), until after 1/Ro 1, N u drops monotonically with increasing 1/Ro (Zhong et al. 2009). Figure  9(a) shows that for span-wise rotating RB convection N u also initially does not change much for 1/Ro 1. After that N u drops monotonically until reaching its minimum at 1/Ro ≈ 10, and then it increases monotonically towards the 2D value for a tworoll state. Similar non-monotonic behavior of N u with the control parameter was also found in sheared RB convection, where N u also non-monotonically depends on the wall shear Reynolds number (Blass et al. 2020). For span-wise rotating RB convection, Re monotonically increases from the 3D value towards the 2D value with increasing rotation rate, as shown in figure 9(b). Note that this behavior is very different from RB convection rotating around the vertical axis, where Re decreases monotonically with increasing 1/Ro (Chong et al. 2017).
We now connect the global transport properties with the flow organization. Figure 10 shows instantaneous temperature fields at the midheight (top row) and boundary layer height close to the bottom plate (bottom row) for different 1/Ro with Ra = 10 7 , P r = 0.71, Γ = 16. We can clearly see the connection between large-scale thermal structure at midheight and boundary layer height for different 1/Ro, which has also been shown in 3D RB convection with no-slip plates (Stevens et al. 2018;Green et al. 2020). With increasing rotation rate, one sees the increasing two-dimensionlization of the flow. For the non-rotation case (1/Ro = 0), figure 10(a) shows qualitatively similar superstructures as the no-slip case (Pandey et al. 2018;Stevens et al. 2018). When 1/Ro increases to 1, a meandering large-scale convection roll state develops, as can be seen from figure 10(b). Interestingly, similar meandering structures have also been observed in many shear-driven flows when the horizontal isotropy is broken, such as plane Couette flow (Lee & Moser 2018), wavy Taylor rolls in Taylor-Couette flow (Andereck et al. 1986), and also in sheared RB convection (Blass et al. 2020). Figure 9(a) shows that this meandering structure at 1/Ro = 1 has still similar N u as in the non-rotation case. For medium rotation rates Figure 9: 3D RB convection with span-wise rotation: (a) N u and (b) Re as function of 1/Ro for Ra = 10 7 , P r = 0.71, Γ = 16 (black circles). For orientation with respect to the Nusselt number, the data for non-rotation (1/Ro = 0, red diamond) and the 2D cases with the same control parameters (Ra = 10 7 , P r = 0.71, Γ = 16) for different roll aspect ratios Γ r (blue squares) are also shown; for these data points the value at the 1/Ro axis has no meaning. The Reynolds number Re for the Γ r = 8/3 (7702.85) and Γ r = 2 (7726.20) states are close to each other and can not be differentiated in the figure. 1/Ro = 3.75 and 10, a two-roll state evolves; interestingly, the cyclonic circulation has a larger size than the anticyclonic one. The smaller N u for these two-roll states is related to the decreased plume emission area. With further increasing rotation rate, the flow increasingly behaves like 2D cases. For the largest 1/Ro = 50 as shown in figure 10(e), a two-roll state with equal size of each roll has developed, with small span-wise variation in temperature.
After we have shown the increasing two-dimensionlization of the flow with increasing rotation rate for span-wise rotating RB convection, next we will study the transition from zonal flow to the convection roll states with increasing aspect ratio Γ , similarly as we have already done for the 2D case. von Hardenberg et al. (2015) studied span-wise rotating RB convection for Ra = 10 7 , P r = 0.71 with fixed Γ = 2π. They observe strong zonal flow that is perpendicular to both rotation vector and gravity vector. Both the cyclonic zonal flow and the anticyclonic one have been obtained using different initial conditions. These two kinds of zonal flow are symmetric for 2D, while they are not in 3D with span-wise rotation, as the Coriolis force depends on the direction of velocity and thus it breaks the symmetry between the two kinds of zonal flow, which have an opposite flow direction. The main difference between the two kinds of zonal flow is that intermittent bursts exist for anticyclonic zonal flow, similar to what is observed in 2D cases with small P r 2, while these bursts are absent for the cyclonic zonal flow (von Hardenberg et al. 2015). We note that, as in von Hardenberg et al. (2015), the dimensionless angular velocity Ω = Ωτ th is used to quantify the rotation velocity, where τ th = H 2 /κ is thermal diffusive time. The dimensionless angular velocity is related to the Rossby number by Ro = √ RaP r/(2Ω ) (Novi et al. 2019). We performed simulations for Ra = 10 7 , P r = 0.71, Γ = 8 and 16 with 1/Ro = 3.75, which corresponds to 2Ω = 10000 in von Hardenberg et al. (2015). Three different initial conditions were used to trigger different possible states, namely: • IC 0 with zero initial velocity, • IC c with cyclonic shear flow u(z) = 2z − 1, v = 0, w = 0, to trigger possible cyclonic zonal flow, and • IC a with anticyclonic shear flow u(z) = 1 − 2z, v = 0, w = 0, to trigger possible anticyclonic zonal flow.
We first report the results for Γ = 8: Figure 11(a) illustrates that for the initial conditions IC 0 , the flow quickly develops into a two-roll state as indicated in figure 11(c).
We note that the cyclonic circulation is again larger than the anticyclonic one. For IC c as initial conditions, the cyclonic zonal flow shown in figure 11(d) can statistically stably exist for more than 3000 free-fall time units. Finally, the initial conditions IC a can trigger an anticyclonic zonal flow with burst phenomenon, which is consistent with the findings of von Hardenberg et al. (2015). However, here this feature only lasts for about 380 freefall time units and then the system undergoes a transition to a two-roll state. We thus conclude that for Γ = 8 the system can again display a bi-stability behavior, in which both zonal flow and convection roll states are statistically stable.
We now come to the case of Γ = 16. From figure 11(b) we conclude that for three different initial conditions, the flow eventually evolves to the very same final state, namely the convection roll state. The cyclonic zonal flow initially seen for short time for the initial conditions IC c quickly undergoes a transition to the convection roll state. The snapshots in figure 11(b − d) reveal similar transition processes as we had already observed in 2D. Again, the final two-roll state has a larger cyclonic circulation. For the final convection state the horizontal scale of the flow reaches the domain size of 16, which is much larger than the typical horizontal scale of superstructures observed in 3D RB convection with no-slip plates (Pandey et al. 2018;Stevens et al. 2018). Such large-scale structures can not be captured in small domains, which is the reason why for small domains only zonal flow states can be realized (von Hardenberg et al. 2015).
To summarize our results from our 3D simulations with free-slip plates and span-wise rotation (with fixed Ra = 10 7 , P r = 0.71, and relative strong rotation 1/Ro = 3.75): We have revealed a similar physical picture as we had done before for 2D RB convection with free-slip plates: (i) For small aspect ratio Γ = 2π, the flow is zonal (von Hardenberg et al. 2015).
(ii) With increasing Γ up to at least Γ = 8, the convection roll state and the zonal flow state coexist in phase space and which one is taken depends on the initial conditions.
(iii) For larger Γ = 16, we always obtain the convection roll state, independently on of what kind of initial conditions were employed.

Concluding remarks
In summary, we have studied 2D RB convection and span-wise rotating RB convection with free-slip plates in horizontally periodic domains using direct numerical simulations. We find that the zonal flow, which was previously observed in a Γ = 2 cell (van der Poel et al. 2014;Goluskin et al. 2014), cannot be sustained and will undergo transitions to convection roll states when the aspect ratio Γ is larger than a critical value, which depends on the Rayleigh number Ra and Prandtl number P r.
We reveal three regimes: (i) For small Γ (typically Γ 1 − 3, depending on Ra and P r), only zonal flow can be observed.
(ii) With increasing Γ , we first find a bi-stable regime in which, depending on the initial conditions, both zonal flow and convection roll states can be stable .
(iii) For even larger-aspect ratio systems, only convection roll states can be sustained. How many convection rolls develop in the convection roll states (in other words, what is the mean aspect ratio Γ r of an individual roll) depends on the initial conditions. For instance, for Ra = 10 8 and P r = 10 the horizontal extent of the stable convection rolls varies between 16/11 and 64 times the height of the convection cell. A convection roll with an as large aspect ratio of Γ r = 64 or more generally already with Γ r 10 can be seen as localized zonal flow.
The heat transfer in the system increases significantly when the horizontal extent of the convection roll is reduced. It is found that the Prandtl number dependence of the Nusselt number N u(P r) has very different trends for large and small Γ r : For large Γ r (like Γ r = 16 for a Γ = 32 cell), the flow behaves like a "localized" zonal flow state, and N u increases with increasing P r, similarly as we found for the zonal flow state. In contrast, for small Γ r , the heat flux into the system is dominated by the plume-impacting regions on the bottom plate, in which the local heat flux is very high and increases with decreasing P r, implying that the global heat flux N u also increases with decreasing P r.
For span-wise rotating 3D RB convection, we find that with increasing rotation rate 1/Ro, both the transport properties (like the Nusselt number N u and the Reynolds number Re) and the flow organization increasingly behave like in the corresponding 2D cases. In particular, just as in 2D, the zonal flow observed in a small periodic cell with Γ = 2π (von Hardenberg et al. 2015), disappears in larger cells with Γ = 16. For intermediate Γ = 8, bistability is observed, again similarly as observed in 2D RB convection.
Finally, an interesting but still open question is the final fate of the aspect ratio dependence of the zonal flow for higher Ra, i.e., is there a finite Ra above which zonal flow exists for all Γ ? Due to the chaotic nature of the flow, mapping out the parameter regime where zonal flow can be found is not easy, especially not for high Ra and large Γ .
From a broader perspective, our study underlines the importance of having large enough aspect ratios in numerical simulations of wall-bounded turbulent flows, even when one employs periodic boundary conditions. We had seen this before in 3D RB convection with no-slip velocity boundary conditions at the plates (Pandey et al. 2018;Stevens et al. 2018;Krug et al. 2020;Green et al. 2020), but apparently this conclusion is much more general. Table 3: Simulation details for all cases shown in figure 6(a). The columns from left to right indicate Ra, P r, Γ , grid resolutions N x × N z , the number of initial rolls n (i) , the number of final convection rolls n, the mean aspect ratio of the convection rolls Γ r = Γ/n, the Nusselt number N u, the Reynolds number Re based on root-mean-square of the global velocity, the horizontal Reynolds number Re x based on root-mean-square of the horizontal velocity, the vertical Reynolds number Re z , the total simulation time t tot , and the time t avg used to average N u and Re.   Table 5: Simulations details for 3D RB convection with span-wise rotation. The corresponding 2D simulations are also included for comparison. The table head is similar to that of table 3 and 4, apart from the inverse Rossby number 1/Ro for the 3D cases. IC 0 means initial conditions with zero velocity and conductive temperature profile with random perturbations. IC c means initial conditions with cyclonic shear flow u(z) = 2z−1, v = 0, w = 0 and conductive temperature profile with random perturbations. Figure 12: (a) Reynolds number ratio Re z /Re x as function of Γ for the zonal flow state and (b) Reynolds number ratio Re z /Re x as function of Γ r for the convection roll state for Ra = 10 8 with different P r. Solid symbols in (b) for Γ = 16 and hollow symbols for Γ = 32. The solid symbols can hardly be seen as they are mostly hidden by the hollow symbols. The data in (b) can be well described by the effective scaling relation Re z /Re x = 0.86Γ −0.68 r . Appendix B. Reynolds number ratio Figure 12(a) shows the Reynolds number ratio Re z /Re x as function of Γ for the zonal flow state. The ratio Re z /Re x increases with increasing P r for P r 10, which has a similar trend as N u discussed before. Figure 12(b) shows Reynolds number ratio Re z /Re x as function of Γ r for convection roll states for Ra = 10 8 with different P r. It is remarkable that Re z /Re x seems to have a universal dependence on Γ r for different P r, despite that N u has a different trend on P r for large and small Γ r . The data can be well represented by the effective scaling relation Re z /Re x = 0.86Γ −0.68 r .