Regime transitions in thermally driven high-Rayleigh number vertical convection

Abstract Thermally driven vertical convection (VC) – the flow in a box heated on one side and cooled on the other side, is investigated using direct numerical simulations with Rayleigh numbers over the wide range of $10^7\le Ra\le 10^{14}$ and a fixed Prandtl number $Pr=10$ in a two-dimensional convection cell with unit aspect ratio. It is found that the dependence of the mean vertical centre temperature gradient $S$ on $Ra$ shows three different regimes: in regime I ($Ra \lesssim 5\times 10^{10}$), $S$ is almost independent of $Ra$; in the newly identified regime II ($5\times 10^{10} \lesssim Ra \lesssim 10^{13}$), $S$ first increases with increasing $Ra$ (regime $\textrm {{II}}_a$), reaches its maximum and then decreases again (regime $\textrm {{II}}_b$); and in regime III ($Ra\gtrsim 10^{13}$), $S$ again becomes only weakly dependent on $Ra$, being slightly smaller than in regime I. The transition from regime I to regime II is related to the onset of unsteady flows arising from the ejection of plumes from the sidewall boundary layers. The maximum of $S$ occurs when these plumes are ejected over approximately half of the area (downstream) of the sidewalls. The onset of regime III is characterized by the appearance of layered structures near the top and bottom horizontal walls. The flow in regime III is characterized by a well-mixed bulk region owing to continuous ejection of plumes over large fractions of the sidewalls, and, as a result of the efficient mixing, the mean temperature gradient in the centre $S$ is smaller than that of regime I. In the three different regimes, significantly different flow organizations are identified: in regime I and regime $\textrm {{II}}_a$, the location of the maximal horizontal velocity is close to the top and bottom walls; however, in regime $\textrm {{II}}_b$ and regime III, banded zonal flow structures develop and the maximal horizontal velocity now is in the bulk region. The different flow organizations in the three regimes are also reflected in the scaling exponents in the effective power law scalings $Nu\sim Ra^\beta$ and $Re\sim Ra^\gamma$. Here, $Nu$ is the Nusselt number and $Re$ is the Reynolds number based on maximal vertical velocity (averaged over vertical direction). In regime I, the fitted scaling exponents ($\beta \approx 0.26$ and $\gamma \approx 0.51$) are in excellent agreement with the theoretical predictions of $\beta =1/4$ and $\gamma =1/2$ for laminar VC (Shishkina, Phys. Rev. E., vol. 93, 2016, 051102). However, in regimes II and III, $\beta$ increases to a value close to 1/3 and $\gamma$ decreases to a value close to 4/9. The stronger $Ra$ dependence of $Nu$ is related to the ejection of plumes and the larger local heat flux at the walls. The mean kinetic dissipation rate also shows different scaling relations with $Ra$ in the different regimes.


Introduction
Thermally driven convective fluid motions are ubiquitous in various geophysical and astrophysical flows, and are important in many industrial applications. Rayleigh-Bénard convection (RBC) (Ahlers, Grossmann & Lohse 2009;Lohse & Xia 2010;Chillà & Schumacher 2012;Xia 2013), where a fluid layer in a box is heated from below and cooled from above, and vertical convection (VC) (Ng et al. 2015;Shishkina 2016;Ng et al. 2017Ng et al. , 2018, where the fluid is confined between two differently heated isothermal vertical walls, have served as two classical model problems to study thermal convection. Vertical convection was also called convection in a differentially heated vertical box in many early papers (Paolucci & Chenoweth 1989;Le Quéré & Behnia 1998). Both RBC and VC can be viewed as extreme cases of the more general so-called tilted convection (Guo et al. 2015;Shishkina & Horn 2016;Wang et al. 2018a,b;Zwirner & Shishkina 2018;Zwirner et al. 2020;Zhang, Ding & Xia 2021), with a tilt angle of 0 • for RBC and 90 • for VC. We focus on VC in this study. Vertical convection finds many applications in engineering, such as thermal insulation using double-pane windows or double walls, horizontal heat transport in water pools with heated/cooled sidewalls, crystal growth procedures, nuclear reactors, ventilation of rooms, and cooling of electronic devices, to name only a few. Vertical convection has also served as a model to study thermally driven atmospheric circulation (Hadley 1735;Lappa 2009) or thermally driven circulation in the ocean, e.g. next to an ice-block (Thorpe, Hutt & Soulsby 1969;Tanny & Tsinober 1988).
The main control parameters in VC are the Rayleigh number Ra ≡ gαL 3 Δ/(νκ) and the Prandtl number Pr ≡ ν/κ. Here, α, ν and κ are the thermal expansion coefficient, the kinematic viscosity and the thermal diffusivity of the convecting fluid, respectively, g is the gravitational acceleration, Δ ≡ T h − T c is the temperature difference between the two side walls, and L is the width of the convection cell. The aspect ratio Γ ≡ H/L is defined as the ratio of height H over width L of the domain. The responses of the system are characterized by the Nusselt number Nu ≡ QL/(kΔ) and the Reynolds number Re ≡ UL/ν, 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 is the characteristic velocity of the flow.
Since the pioneering work of Batchelor (Batchelor 1954), who first addressed the case of steady-state heat transfer across double-glazed windows, VC has drawn significant attention especially in the 1980s and 1990s, and most of these studies used experiments or a two-dimensional (2-D) direct numerical simulation (DNS) in a square domain with unit aspect ratio. For relatively low Ra (e.g. Ra < 10 3 ), the flow is weak and heat is transferred mainly by thermal conduction. With increasing Ra, typical stratified flow structures appear in the bulk region (de Vahl Davis & Jones 1983), while the flow remains steady. With a further increase in Ra, the flow becomes unsteady with periodical/quasi-periodical or chaotic motions (Paolucci & Chenoweth 1989;Le Quéré & Behnia 1998), and eventually becomes turbulent when Ra is sufficiently high (Paolucci 1990).
The onset of unsteadiness has been well explored in the past (Chenoweth & Paolucci 1986;Paolucci & Chenoweth 1989;Janssen & Henkes 1995;Le Quéré & Behnia 1998). Paolucci & Chenoweth (1989) investigated the influence of the aspect ratio Γ on the onset of unsteadiness for 2-D VC with Pr = 0.71. They found that for Γ 3, the first transition from the steady state arises from an instability of the sidewall boundary layers, while for smaller aspect ratios 0.5 ≤ Γ 3, it arises from internal waves near the departing corners. Such oscillatory instability arising from internal waves was first reported by Chenoweth & Paolucci (1986). Paolucci & Chenoweth (1989) also found that for Γ = 1, the critical Rayleigh number Ra c for the onset of unsteadiness lies between 1.8 × 10 8 and 2 × 10 8 . Later work, with Pr = 0.71 and Γ = 1 by Le Quéré & Behnia (1998), also showed that the internal gravity waves play an important role in the time-dependent dynamics of the solutions, and 1.81 × 10 8 ≤ Ra c ≤ 1.83 × 10 8 was determined to be the range of the critical Rayleigh number. Janssen & Henkes (1995) studied the influence of Pr on the instability mechanisms for Γ = 1, and observed that for 0.25 ≤ Pr ≤ 2, the transition occurs through periodic and quasi-periodic flow regimes. One bifurcation is related to an instability occurring in a jet-like fluid layer exiting from the corners of the cavity where the vertical boundary layers are turned horizontal. Such jet-like flow structures are responsible for the generation of internal gravity waves (Chenoweth & Paolucci 1986;Paolucci & Chenoweth 1989). The other bifurcation occurs in the boundary layers along the vertical walls. Both of these instabilities are mainly shear-driven. For 2.5 ≤ Pr ≤ 7, Janssen & Henkes (1995) found an 'immediate' (i.e. sharp) transition from the steady to the chaotic flow regime, without intermediate regimes. This transition is also caused by boundary layer instabilities. They also showed that Ra c significantly increases with increasing Pr, e.g. for Pr = 4, the flow can still be steady with Ra = 2.5 × 10 10 . However, owing to the computation limit, unsteady motions for the large-Pr cases have largely remained unexplored in the past.
Additionally, the flow structures for VC have been examined in detail. A typical flow feature for VC is the stably-stratified bulk region (de Vahl Davis & Jones 1983;Ravi, Henkes & Hoogendoorn 1994;Trias et al. 2007;Sebilleau et al. 2018;Chong et al. 2020). Such stratification can be quantified by the time-averaged non-dimensional temperature gradient at the centre, namely (1.1) Here t denotes a time average. Gill (1966) derived asymptotic solutions for high Pr, and predicted S = 0.42 as Ra → ∞, while an accurate solution of the same system by Blythe, Daniels & Simpkins (1983) predicted a value of 0.52. Later DNS results for Ra = 10 8 and Pr = 70 yielded S = 0.52 (Ravi et al. 1994), which is in excellent agreement with the theoretical prediction by Blythe et al. (1983). However, for small Pr, the structure of the core and the vertical boundary layer are no longer similar to those predicted by the asymptotic solutions which are valid for large Pr (Blythe et al. 1983). Unfortunately, there exists no such asymptotic theory for finite Pr. Only Graebel (1981) has presented some approximate solutions, in which some terms have been neglected in the equations. For Pr = 0.71, his prediction yielded S = 0.49, which is considerably smaller than the value S ≈ 1 from DNS (Ravi et al. 1994;Trias et al. 2007). It was concluded that S is independent of Ra for Ra ≤ 10 10 (Paolucci 1990); however, it is evident that the dependence of S on Ra and Pr, especially for those with high Ra > 10 10 and low Pr < 0.71, is still poorly understood.
A key question in the study of thermal convection is: How do Nu and Re depend on Ra and Pr? This question has been extensively addressed in RBC over the past years (Ahlers et al. 2009). For RBC, the mean kinetic dissipation rate ( u ) and thermal dissipation rate ( θ ) obey exact global balances, which feature Ra, Nu and Pr (Shraiman & Siggia 1990). For this problem, in a series of papers, Grossmann & Lohse (2000, 2001, 2002, 2004) developed a unifying theory to account for Nu (Ra, Pr) and Re(Ra, Pr) over wide parameter ranges. The central idea of the theory is a decomposition of u and θ into their boundary layer and bulk contributions. The theory has been well confirmed through various experiments and numerical simulations (Stevens et al. 2013). This theory has also been applied to horizontal convection (Shishkina, Grossmann & Lohse 2016;Shishkina & Wagner 2016) and internally heated convection (Wang, Shishkina & Lohse 2020b). However, in VC, the exact relation for u does not hold, which impedes the applicability of the unifying theory to the scalings in VC (Ng et al. 2015).
Most of the simulations for VC were conducted for Ra 10 10 . The high-Ra simulations become stiff owing to a decrease in the boundary-layer thicknesses with increasing Ra. As a result, little is known about what will happen at Ra much larger than 10 10 . In this study, we attempt to fill this gap in knowledge by performing DNS up to Ra = 10 14 . The price we have to pay is that for such large Ra, we are restricted to 2-D.
The main questions we want to address in this study are as follows.
(i) Is the conclusion that S is independent of Ra for Ra ≤ 10 10 (Paolucci 1990) still valid for Ra much larger than 10 10 ? (ii) How does the global flow organization (mean temperature and velocity profiles) change with increasing Ra up to 10 14 ? (iii) How robust are the laminar scaling relations Nu ∼ Ra 1/4 and Re ∼ Ra 1/2 (Shishkina 2016) for higher Ra? Will new scaling relations appear for Ra much larger than 10 10 ?
We find that S is not independent of Ra over the studied parameter range at all. Instead, we find that apart from the small-Ra regime (now called regime I), where S only weakly depends on Ra (Paolucci 1990), there are additional regimes for Ra 5 × 10 10 with different scaling relations. In regime II (5 × 10 10 Ra 10 13 ), with increasing Ra, S first increases (regime II a ) to its maximum and then decreases (regime II b ) again. In regime III (Ra 10 13 ), S again becomes weakly dependent on Ra, with a smaller value than that of regime I. Furthermore, we find that the laminar power law exponents β = 1/4 and γ = 1/2 undergo sharp transitions to β ≈ 1/3 and γ ≈ 4/9 when Ra 5 × 10 10 , i.e. at the transition from regime I to regime II.
The rest of the paper is organized as follows. Section 2 describes the governing equations and numerical methods. The different flow organizations in the different regimes are studied in § 3. In § 4, we discuss the transition of the scaling relations for heat and momentum transport between the different regimes. Finally, § 5 contains a summary and an outlook.

Numerical procedures
A sketch of 2-D VC is shown in figure 1. The top and bottom walls are insulated. The left wall is heated with temperature T h , while the right wall is cooled with temperature T c . No-slip and no-penetration velocity boundary conditions are used at all the walls. The aspect ratio Γ ≡ H/L is fixed to 1. The dimensionless governing equations are the incompressible Navier-Stokes equations with an Oberbeck-Boussinesq approximation: Here e z is the unit vector pointing in the direction opposite to gravity. The dimensionless velocity, temperature and pressure are represented by u ≡ (u, w), θ and p, respectively. For non-dimensionalization, we use the width of the convection cell L and the free-fall velocity U = (gαΔL) 1/2 . Temperature is non-dimensionalized as θ = (T − T c )/Δ. The governing equations were solved using the second-order staggered finite-difference code AFiD (Verzicco & Orlandi 1996;van der Poel et al. 2015). The code has already been extensively used to study RBC (Wang et al. 2020a,c;Liu et al. 2021) and internally heated convection (Wang et al. 2020b). Direct numerical simulation was performed for 10 7 ≤ Ra ≤ 10 14 with a fixed Pr = 10. Stretched grids were used to resolve the thin boundary layers and adequate resolutions were ensured to resolve the small scales of turbulence (Shishkina et al. 2010). Grids with up to 8192 × 8192 nodes were used for the highest Ra = 10 14 . We performed careful grid independence checks for several high-Ra cases. It was found that the difference of Nu and Re for the different grids were always smaller than 1 % and 2 %, respectively. Details on the simulations are provided in table 2 in the appendix.

Global flow fields
We first focus on the change of global flow organizations with increasing Ra. Figure 2 shows instantaneous temperature, horizontal velocity (u) and vertical velocity (w) fields 917 A6-5 for different Ra. For the considered Pr = 10, we find that the flow is still steady for Ra = 5 × 10 10 , as shown in figure 2(a-c), which is consistent with the finding that the critical Rayleigh number Ra c for the onset of unsteadiness increases with increasing Pr and that the flow is indeed still steady for Ra = 2.5 × 10 10 with Pr = 4 (Janssen & Henkes 1995). This is in sharp contrast with RBC, where the flow is already turbulent for such high Ra with Pr = 10 ( Wang et al. 2020c). The flow is stably stratified in the bulk region, as shown in figure 2(a). The large horizontal velocity regions mainly concentrate near the top and bottom walls (figure 2b), while the strong vertical motion mainly occurs near the two sidewalls (figure 2c). Such flow structures are typical for steady VC with large Pr (Ravi et al. 1994). However, with a minor increase of Ra from Ra = 5 × 10 10 to Ra = 6 × 10 10 , the flow becomes instantaneously chaotic, as shown in figure 2(d-f ). This finding is consistent with the previous result that for Pr ≥ 2.5, there is an immediate transition from the steady to the chaotic flow regime without intermediate regimes (Janssen & Henkes 1995). This transition is caused by boundary layer instabilities, which are reflected in the plume ejections in the downstream of the boundary layers (figure 2d). The strong horizontal/vertical fluid motions still concentrate near the horizontal/vertical walls, as indicated in figures 2(e) and 2( f ). However, there are already some chaotic features appearing in the bulk, which suggest a change of the bulk properties.
When Ra is further increased to 6 × 10 11 (figure 2g-i), further evident changes of the global flow organization appear as follows. (i) The hot plumes mainly eject over the upper half of the hot sidewall, and enter the upper half of the bulk region. This makes the hot upper bulk region more isothermal than in the smaller-Ra cases. Similar processes happen for the cold plumes and the lower cold bulk region. Therefore, figure 2(g) clearly shows a larger centre temperature gradient than those in figures 2(a) and 2(d). (ii) The strong horizontal motions now not only occur near the horizontal walls, but also in the bulk region (figure 2h), and alternating rightward and leftward 'zonal flow' structures appear.
For the highest Ra = 10 14 (figure 2j-l), the thermal driving is so strong that hot plumes are now ejected over large fractions of the left vertical wall (0.2 z/L 1). The plumes are transported into the bulk region by the zonal flow structures shown in figure 2(k). This process causes efficient mixing in the bulk, which then leads to a smaller centre

Mean profiles for temperature and horizontal velocity
We have seen that the global flow organization evidently changes with increasing Ra. In this subsection, we quantify these changes by looking at the mean profiles for the temperature and for the horizontal velocity. Figure 3(a) clearly shows the change in the temperature profiles at x/L = 0.5 with increasing Ra, which is consistent with the temperature fields presented in figure 2. The change of the bulk temperature profile shape can be quantified by the time-averaged non-dimensional vertical temperature gradient in the cell centre, i.e. S ≡ (L/Δ)(∂T/∂z) c t (Paolucci 1990;Ravi et al. 1994). This quantity is plotted in figure 3(b), where one can observe three different regimes. In the well-explored regime I (Ra 5 × 10 10 ), S weakly depends on Ra, with a value S ≈ 0.5, close to that of S = 0.52 for Ra = 10 8 with Pr = 70 reported in Ravi et al. (1994). However, in regime II (5 × 10 10 Ra 10 13 ), S has a non-monotonic dependence on Ra: it first increases with increasing Ra, reaches its maximum at Ra = 6 × 10 11 , and then decreases again. Regime II is further divided into regime II a , where S increases with increasing Ra, and regime II b , where S decreases with increasing Ra. The onset of regime II coincides with the onset of unsteadiness, which shows that plume emissions play an important role in altering the bulk properties. The maximum of S occurs when approximately half of sidewall areas (downstream) feature plume emissions, as shown in figure 2(g). Finally, in regime III, S again becomes weakly dependent on Ra, while it has a smaller value than that of regime I. The small value of S in regime III arises from the well-mixed bulk region, as can be seen in figure 2( j). Figure 3(c) shows the change of the horizontal velocity profiles with increasing Ra. In regime I, the strong horizontal fluid motions concentrate near the top and bottom walls. In contrast, in regime III, the largest horizontal velocity appears in the bulk region, and alternating rightward and leftward fluid motions, i.e. zonal flows, are observed even after time averaging, which is consistent with the instantaneous horizontal velocity field shown in figure 2(k). Another prominent flow feature of regime III is that the horizontal velocity near the top and bottom walls is close to 0. This means that the 'layered structure' near the top and bottom walls, as indicated in the temperature field in figure 2( j), is actually nearly a 'dead zone' with weak fluid motions. Thus the appearance of this nearly 'dead' layered structure indicates the the onset of regime III. Regime II serves to connect regime I and regime III: in regime II a , the strongest horizontal motion still takes place near the top and bottom walls, see, e.g. the horizontal velocity profile for Ra = 2 × 10 11 in figure 3(c). In contrast, in regime II b , the strongest horizontal fluid motions appear in the bulk, see, e.g. the strong zonal flow motions in the central region for Ra = 10 12 , as shown in figure 3(c). We remark that the zonal flow has been found in many geo-and astrophysical flows (Yano, Talagrand & Drossart 2003;Heimpel, Aurnou & Wicht 2005;Nadiga 2006), and it has also been extensively studied in RBC (Goluskin et al. 2014;Wang et al. 2020a;Zhang et al. 2020;Reiter et al. 2021). It is remarkable and interesting to also observe zonal flows in the high-Ra VC system. This system thus provides another model to study the physics of the zonal flow.

Mean profiles for temperature and vertical velocity
We now consider the mean vertical velocity and temperature profiles in the longitudinal (x) direction. Figure 4(a) shows the mean temperature profiles in the longitudinal (x) direction at mid-height z/L = 0.5 for different Ra. It is seen that for Ra values that are not too high, the temperature does not monotonically drop from θ = 1 at the hot wall to θ = 0.5 in the core. Instead, an undershoot phenomenon is observed. This phenomenon arises from the stable stratification in the bulk (Ravi et al. 1994) and can also be observed in the similarity solutions of the boundary layer equations for natural convection over a vertical hot wall in a stably stratified environment (Henkes & Hoogendoorn 1989). However, for Ra ≥ 10 12 , we find that the overshoot phenomenon disappears. This arises from the continuous emissions of hot plumes at mid-height z/L = 0.5 and beyond, as then the hot fluid directly touches the well-mixed bulk flow with small stratification. Figure 4(b) shows vertical velocity profiles in the longitudinal direction, again at mid-height z/L = 0.5. With increasing Ra, the boundary layer becomes thinner and the peak vertical velocity becomes smaller. This finding reflects the different flow organizations in the different regimes: the emitted plumes in regimes II and III weaken the overall vertical fluid motions, as compared with the steady flow organization in regime I. This is also reflected in the Re ∼ Ra γ scaling, as will be discussed below.

Global heat and momentum transport and dissipation rates
Next, we focus on the global heat (Nu) and momentum (Re) transport. Here, we use the wind-based Reynolds number Re with the characteristic velocity which is the same definition as that in Shishkina (2016), and the root-mean-square Reynolds number Re rms with the characteristic velocity where V,t indicates volume and time averaging. Figures 5(a) and 5(b) show that in regime I, the obtained effective power law scaling relations agree remarkably well with the theoretical prediction made for laminar VC (Shishkina 2016), namely, Nu ∼ Ra 1/4 and Re ∼ Ra 1/2 . The fitted scaling relations are provided in table 1. It is also seen that a slightly faster growth of Nu with Ra is obtained for Ra ≤ 10 9 . A similar increase of the scaling exponent for small Ra has also been found previously in both confined (Shishkina 2016;Wang, Zhang & Guo 2017;Wang et al. 2019) and double periodic VC (Ng et al. 2015). However, when Ra ≥ 5 × 10 10 , in regime II and regime III, evidently different scaling relations are observed. The fitted power law scaling relations (see theory for RBC (Grossmann & Lohse 2000). Such observation of scaling transitions further demonstrates that there are no pure scaling laws in thermal convection. This has already been seen in RBC (Grossmann & Lohse 2000;Ahlers et al. 2009), horizontal convection Shishkina & Wagner 2016;Reiter & Shishkina 2020) and internally heated convection (Wang et al. 2020b), and apparently crossovers between different scaling regimes also occur here. However, the sharpness of the scaling transition from β = 1/4 to β = 1/3 observed here is quite different from the smooth transition seen in RBC. Indeed, in RBC, the transition from Nu ∼ Ra 1/4 to Nu ∼ Ra 1/3 is very smooth, spread over more than two orders of magnitude in Ra (Grossmann & Lohse 2000), and the linear combination of the 1/4 and 1/3 power laws even mimics an effective 2/7 scaling exponent (Castaing et al. 1989) over many orders of magnitude in Ra. Figure 5(c) shows that also Re rms behaves differently in different regimes. The fitted scaling relation Re rms ∼ Ra 0.37 in regime I is the same as that found for 2-D VC with Pr = 0.71 (Wang et al. 2019), which suggests that in vertical convection, different Reynolds numbers have different scaling relations with Ra. In regime II, the normalized Reynolds number Re rms /Ra 1/2 depends non-monotonically on Ra, and shows a pronounced local minimum. The fitted scaling exponent 0.31 in regime III is again smaller than that in regime I.
It is interesting to note that, though S shows a clear transition between regime II and regime III, this transition is not seen in the effective power law scaling relations Nu ∼ Ra β and Re ∼ Ra γ . Our interpretation of this noteworthy finding is as follows. The transition of the flow organization from regime I to regime II is sharp in view of the sudden appearance of plume emissions from the sidewall thermal boundary layers, and in view of the emergence of the local minimum of the Nusselt number distribution on the sidewall. However, the transition of flow organization from regime II to regime III is more continuous. The flows in these two regimes are characterized by the alternating rightward and leftward fluid motions, i.e. zonal flows, in the bulk, and they all have plume emissions from the sidewall boundary layers and a local minimum of the Nusselt number distribution on the sidewall. The only prominent difference is the appearance of layered structures near the top and bottom plates in regime III. As this layered structure only concentrates in a small region near the top and bottom plates, the global effective scaling exponents of Nu and Re do not seem to be sensitive to the different flow organizations in regimes II and III.
To better understand the sudden change of the global heat transport properties at the transition to regime II, we now consider the wall heat flux, which is denoted by the local Nusselt number at the wall Nu(z)| x=0,1 = ∂ θ t /∂x| x=0,1 . Figure 6(a) displays Nu(z)| x=0 at the left wall, while Nu(z)| x=1 at the right wall is not shown owing to the inherent symmetry of the system. For Ra ≤ 5 × 10 10 , the local Nu(z)| x=0 generally decreases monotonically with increasing heights z. The large local Nu(z)| x=0 for small heights z is attributed to the fact that the hot fluid there is in direct contact with the cold fluid, which leads to large temperature gradients. In contrast, for Ra > 5 × 10 10 , a local minimum 917 A6-11 and a local maximum in Nu(z)| x=0 are identified, the heights of which are denoted as z t1 and z t2 . It is clearly seen that Nu(z)| x=0 after the first transition point z t1 increases compared with the steady cases at Ra ≤ 5 × 10 10 . This is because the emissions of the plumes lead to more efficient shear-driven mixing, and therefore larger local Nu(z)| x=0 . Thus, the overall heat transport also increases in regime II and later regime III (figure 5a) owing the ejections of plumes, and the change of the scaling is also related to the change of the boundary layer properties. We have shown that the two transition points z t1 and z t2 roughly correspond to the locations where plumes begin to be ejected. Figure 6(b) shows, as expected, that both z t1 and z t2 decrease with increasing Ra, which suggests that the locations where hot plumes begin to be ejected move downwards with increasing Ra. At Ra = 6 × 10 11 , where the centre temperature gradient S achieves its maximum, it can be seen that the mid-height z/L = 0.5 lies between the two transition points, further demonstrating that the maximum of S is achieved once plumes are ejected over approximately half of the area (downstream) of the sidewalls, as seen in the temperature field in figure 2(g).
Finally, we discuss the thermal and kinetic dissipation rates. In RBC, the following exact relations hold (Shraiman & Siggia 1990). The average V is over the whole volume and over time. In VC, the relation (4.4) still holds, however, the relation (4.3) does not hold any longer. Following Ng et al. (2015) and Reiter & Shishkina (2020), we decompose the dissipation rates into their mean and fluctuating parts as Figure 7(a) shows that the relation (4.4) is fulfilled in the DNS. It is also seen that the contribution from the mean field decreases with increasing Ra, which suggests that with increasing Ra, turbulent fluctuations play an increasingly more important role on the mixing process.
The kinetic dissipation rate is displayed in figure 7(b). One can see that the values u /[L −4 ν 3 (Nu − 1)RaPr −2 ] are always smaller than the corresponding value, as occurs in RBC. This was already seen in 3-D VC (Shishkina 2016). For the steady VC with Ra ≤ 5 × 10 10 , the normalized kinetic dissipation rate only weakly depends on Ra, as has also been found in 3-D VC (Shishkina 2016). However, in regimes II and III, it is observed that the normalized kinetic dissipation rate decreases much faster than that in regime I. This can be related to the fact that Nu increases faster in regimes II and III than in regime I. It is also seen that the contribution from turbulent fluctuations is small, similar as in horizontal convection (Reiter & Shishkina 2020).

Conclusions
In conclusion, we have studied vertical convection by direct numerical simulations over seven orders of magnitude of Rayleigh numbers, i.e. 10 7 ≤ Ra ≤ 10 14 , for a fixed Prandtl number Pr = 10 in a two-dimensional convection cell with unit aspect ratio. The main conclusions, which correspond to the answers of the questions posed in the introduction, are summarized as follows.
(i) The dependence of the non-dimensional mean vertical temperature gradient at the cell centre S on Ra shows three different regimes. In regime I (Ra 5 × 10 10 ), S is almost independent of Ra, which is consistent with previous work (Paolucci 1990). However, in the newly identified regime II (5 × 10 10 Ra 10 13 ), S first increases with increasing Ra, reaches its maximum and then decreases again. In regime III (Ra 10 13 ), S again becomes weakly dependent on Ra, with a smaller value than that of regime I. The transition from regime I to regime II coincides with the onset of unsteady fluid motions. The maximum of S occurs when plumes are ejected over approximately half of the area of the sidewall, namely, in the downstream region. The flow in regime III is characterized by a well-mixed bulk region owing to continuous ejection of plumes over large fractions of the sidewalls. Thus S is smaller than that of regime I. (ii) The flow organizations in the three different regimes are quite different from each other. In regime I, the maximal horizontal velocity concentrates near the top and bottom walls. However, the flow gives way to alternating rightward and leftward zonal flows in regime III, where the maximal horizontal velocity appears in the bulk region. Another characteristic feature of the flow in regime III are the 'layered' structures near the top and bottom walls, where the fluid motions are weak. Regime II serves to connect regime I and regime III: in regime II a , the maximal velocity still occurs near the top and bottom walls. In contrast, in regime II b , the zonal flow structures become more pronounced, and the maximal horizontal velocity is found in the bulk region. (iii) Transitions in the scaling relations Nu ∼ Ra β and Re ∼ Ra γ are found. In regime I, the fitted scaling exponents (β ≈ 0.26 and γ ≈ 0.51) are in excellent agreement with the theoretical prediction of β = 1/4 and γ = 1/2 for the laminar VC (Shishkina 2016). However, β increases to a value close to 1/3 and γ decreases to a value close to 4/9 in regimes II and III. The increased heat transport Nu in regimes II and III is related to the ejection of plumes and larger local heat flux at the sidewalls. The mean kinetic dissipation rate also shows different scalings in the different regimes.
We note that the present study only focuses on Pr = 10. Further studies, both numerical simulations and experiments, are needed to address the influence of the Prandtl number Pr and the aspect ratio Γ on the regime transitions in the high-Ra vertical convection. The reported scaling relations for Nu ∼ Ra β and Re ∼ Ra γ and the observed transitions are however already an important ingredient to consider to develop a unifying scaling theory over a broad range of control parameters for vertical convection, to finally arrive at the full dependences Nu(Ra, Pr) and Re(Ra, Pr) and their theoretical understanding. , the Reynolds number based on root-mean-square velocity Re rms , the time t avg used to average Nu and Re. The aspect ratio is fixed to 1 for all the cases. 's' means that the flow is steady. Cases indicated in blue and italic are used for grid independence checks. We note that the difference of Nu for two different grids is always smaller than 1 %, and the difference of Re for the different grids is always smaller than 2 %.