Classical and symmetrical horizontal convection: detaching plumes and oscillations

Classical and symmetrical horizontal convection is studied by means of direct numerical simulations for Rayleigh numbers $Ra$ up to 3 × 1012 and Prandtl numbers $Pr=0.1$, 1 and 10. For both set-ups, a very good agreement in global quantities with respect to heat and momentum transport is attained. Similar to Shishkina & Wagner (Phys. Rev. Lett., vol. 116, 2016, 024302), we find Nusselt number $Nu$ versus $Ra$ scaling transitions in a region $10^{8}\leqslant Ra\leqslant 10^{11}$. Above a critical $Ra$, the flow undergoes either a steady–oscillatory transition (small $Pr$) or a transition from steady state to a transient state with detaching plumes (large $Pr$). The onset of the oscillations takes place at $Ra\,Pr^{-1}\approx 5\times 10^{9}$ and the onset of detaching plumes at $Ra\,Pr^{5/4}\approx 9\times 10^{10}$. These onsets coincide with the onsets of scaling transitions.


Introduction
In a horizontal convection (HC) system, heating and cooling take place over a single horizontal surface of a fluid layer. Sandstrom (1908) argued that, due to the absence of a pressure gradient, a closed circulation cannot be maintained in such systems. However, six decades later, Rossby (1965) demonstrated in his experiments that HC alone, independent of any other sources, is able to create a circulation of a fluid and therefore a net convective buoyancy flux. Over the past decades, Rossby's set-up has become a popular paradigm case to study this important type of natural convective system (Hughes & Griffiths 2008;Griffiths, Hughes & Gayen 2013), which is relevant in geophysical flows like the meridional overturning circulation in the ocean (Munk 1966;Killworth 1983;Scott, Marotzke & Adcroft 2001;Cushman-Roisin & Beckers 2011;Scotti & White 2011), in astrophysical flows (Spiegel 1971) and in engineering applications (Gramberg, Howell & Ockendon 2007). Investigations of HC systems are also necessary for understanding the effect of polar amplification on the ocean circulation (Holland & Bitz 2003), i.e. a phenomenon of global warming that decreases the temperature contrast between the poles and mid-latitudes. In any convective system, a naturally arising question is: How do the global heat transport (Nusselt number Nu) and momentum transport (Reynolds number Re) depend on the main input parameters (Rayleigh number Ra and Prandtl number Pr)? While considering a laminar boundary layer (BL) and balancing buoyancy and viscous dissipation terms inside the BL, Rossby (1965) proposed a relation Nu ∼ Ra 1/5 Pr 0 . The existence of the ∼Ra 1/5 regime was supported by various numerical and experimental studies (e.g. Mullarney, Griffiths & Hughes 2004;Gayen, Griffiths & Hughes 2014), but the predicted Pr invariance of Nu does not hold for small Pr (e.g. Shishkina & Wagner 2016). By considering the dynamics to be driven by a turbulent endwall plume, Hughes et al. (2007) proposed the scaling Nu ∼ Ra 1/5 Pr 1/5 , but, as shown in Shishkina & Wagner (2016), the proposed Pr scaling is too strong and is not supported by direct numerical simulations (DNS). Whereas the Rossby model is based solely on the BL scalings, the model by Shishkina, Grossmann & Lohse (2016) (SGL model), which is an extension of the Grossmann & Lohse (2000 theory to HC, is able to account for laminar regimes as well as for regimes where the mixing is governed by turbulent processes. In particular, the SGL model suggests Nu ∼ Pr 0 Ra 1/4 for large-Pr and Nu ∼ Pr 1/10 Ra 1/5 for small-Pr laminar flows, which was supported by several numerical studies (Shishkina & Wagner 2016;Ramme & Hansen 2019).
However, verification of the other regimes needs further investigations. For this, high-Ra DNS or experiments are needed, which turn out to be challenging tasks. On the one hand, the very slow diffusion in the system is a critical problem for the DNS of (almost) steady flows. On the other hand, in experiments, unwanted heat losses through the vertical walls can affect the scaling results significantly (Ahlers 2000). Therefore, in this work we focus on two set-ups: classical horizontal convection (CHC) and symmetrical horizontal convection (SHC), which can be more suitable for future experiments. Here we report three-dimensional (3-D) DNS results for Ra 3 × 10 12 and Pr = 0.1, 1 and 10.

Theoretical background
We consider a fluid layer that is confined in a rectangular box and heated and cooled locally from the bottom. In the CHC set-up, heating and cooling are applied to the opposite bottom ends (figure 1a); while in the SHC set-up, the bottom is cooled at the ends and heated in the central part (figure 1b). The advantage of the SHC set-up over the HC one is the absence of the vertical endwall attached to the heated plate, which is difficult to isolate thermally in experiments. For both set-ups, we conducted extensive 3-D DNS over the range of parameters shown in figure 1(c), where the DNS data for CHC and Ra < 3 × 10 11 are taken from Shishkina (2017).
where H is the cell height, ν the kinematic viscosity, α the isobaric thermal expansion coefficient, κ the thermal diffusivity and g the acceleration due to gravity. In the CHC configuration, the temperature boundary conditions (BCs) at the bottom are θ = 0.5 for 0 x 0.1 and θ = −0.5 for 0.9 x 1. The other walls are adiabatic, ∂θ /∂n = 0, where n is the wall-normal vector. The velocity BCs are no-slip everywhere. In the SHC set-up, the small vertical endwall near the heated plate is removed and the whole cell is extended by reflection of the cell with respect to the removed endwall. The used finite-volume code is Goldfish (Kooij et al. 2018). A list of all simulations, their spatial resolutions and averaging times are included in the supplementary material available at https://doi.org/10.1017/jfm.2020.211.
The SGL model proposes different scaling regimes based on an assumption that, in HC, the globally averaged kinetic ( u ) and thermal ( θ ) dissipation rates, (2.1) are determined by either the BLs (laminar flows) or the bulk (turbulent flows). Here · V denotes the time and volume average and · z=0 the time and area average at z = 0. All this leads to different scaling regimes of Nu and Re versus Ra and Pr (see table 1).

Global heat and momentum transport
We start our analysis with the Ra dependences of Nu and Re, using the definitions where |∂ z θ c | z=0 is the magnitude of the heat flux considering a pure conductive system subjected to the same BCs (here 1 2 |∂ z θ c | z=0 ≈ 1.12) and Re is based on the total kinetic energy. The results are presented in figure 2. First we observe that Nu and Re in CHC (solid black) and SHC (open colour) match remarkably well, with nearly equal absolute values over the whole parameter range. Therefore, both set-ups can be used for the investigation of the global heat and momentum transport in HC. However, there exist differences in the flow structures, which will be discussed in § 3.2.
The DNS reveal a rather complex scaling dependence, with multiple transitions. Starting from left to right in figure 2(a) we find the following. First, Nu ∼ Ra 1/4 for lower Ra, which corresponds to regime I * l in the SGL model, previously supported by Shishkina & Wagner (2016) and Ramme & Hansen (2019). As Ra increases, all three sets of data for different Pr show a rather sharp transition to a scaling Nu ∼ Ra 1/5 . Note that the critical Ra, where this transition occurs, increases with increasing Pr. The ∼Ra 1/5 scaling seems to persist up to our highest Ra for Pr = 0.1. Regime II l of the SGL model was not observed in our DNS, because Pr = 0.1 is still too large for this regime (Passaggia, Scotti & White 2018). For Pr = 1 and 10, the curves rise again at higher Ra, leading to a scaling exponent of approximately 0.24 and 0.23. Figure  2(b-d) shows Re ∼ Ra 2/5 for small Ra (regime I l ) and Re ∼ Ra 1/3 for larger Ra and no I * l regime at low Pr. However, when Re is based on U 2 + , where · + denotes average in time and over the volume exclusively above the heated plate(s), we find scaling transitions consistent with the Nu-Ra transitions of figure 2(a). In general, the Re scaling is sensitive to its spatial averaging domain. This displays the inhomogeneous nature of HC flows. To explain the scaling transitions, we will have a further closer look at the flow topology and its changes and relate them to the transitions in the scaling relations.

Dynamics: plumes and oscillations
In general, the HC dynamics are rich in flow structures and instability transitions. Paparella & Young (2002) observed that HC flows become unsteady with growing Ra, while higher-Pr flows are stable over a broader range of Ra. Chiu-Webster, Hinch & Lister (2008) and Ramme & Hansen (2019) noticed the existence of time-dependent flows for highly viscous flows. Gayen et al. (2014) showed, for Pr = 5 and varying Ra, that the flow goes through a sequence of stability transitions, starting with the growth of plumes in the BL, followed by convective rolls at higher Ra, and finally show fully 3-D turbulence within a region above the hot BL at Ra ≈ 5 × 10 11 . The linear stability analysis of Passaggia, Scotti & White (2017) for Pr = 1 supports these findings and suggests that there exists a competition between 3-D rolls around a streamwise axis and two-dimensional (2-D) Rayleigh-Taylor (RT) instability. The former seem to be dominant in wider cells, as found by Sukhanovsky et al. (2012), whereas the latter seems to be most relevant for no-slip BC in narrow cells. Sheard & King (2011) found that the onset to unsteady flows is independent of vertical confinement for 0.16 Γ 2, and Passaggia et al. (2018) observed the maximum growth for the 2-D plume instabilities at Γ = 6. HC instability was studied also for the cases of a BL synthetic jet (Leigh, Tsai & Sheard 2016), 2-D HC ) and different temperature BCs (Tsai et al. 2020).
In our DNS we found the existence of 2-D RT instabilities, which manifest themselves as sheared plumes that arise above the heated plate and which travel towards the endwall (CHC) or the centre (SHC), as shown in figure 3(a). However, for small Pr and especially in SHC flows, we found a different time-dependent behaviour prior to plumes emerging, which is an oscillatory instability, that breaks symmetry in SHC (3b). These plume-and oscillatory-induced transitions to the unsteady state are explained below.

Plumes
In terms of time scales, detached plumes can occur if the time scale of the development of RT instabilities T RT is shorter than the advection time scale T wind of the large-scale wind: T RT /T wind < C p , for a certain constant C p . The e-folding time scale of RT instabilities (a characteristic time scale for RT instabilities to grow by the factor 1/e) equals T RT ∼ ν 1/3 /(αg∆) 2/3 (Chandrasekhar 1981) and the time scale of the wind velocity equals T wind ∼ L 2 /(Re ν), which leads to an estimate T wind /T RT = Ra 2/3 /(Re Pr 2/3 ). (3.1) As the plumes detaching regime is anticipated for large Pr, we make use of the scaling relation Re ∼ Pr −1 Ra 1/2 of the regime I * l (see table 1), which gives T wind /T RT ∼ Pr 1/3 Ra 1/6 . Thus, a certain critical value of Pr 1/3 Ra 1/6 , or an equivalent critical value of Pr 2 Ra, determines the onset of the detached plumes. Note that the absolute value of the constant C p can be determined from simulations or experimental data. This relation shows that, for low Pr, the critical Ra increases. Physically explained, the larger wind speed of low-Pr flows advects growing plumes faster to the endwall (CHC) or the centre (SHC) before they become distinguishable from the thermal BL.
The solid red curves for Pr 2 Ra ≈ 10 11 in figure 4(a) and 4(b) give a rough estimate of the Pr and Ra dependence of the onset of the plume-dominated regime. However, using Re scaling relations from the DNS instead of the SGL model, namely Re ∼ Ra 2/5 (figure 2b) and Re ∼ Pr −1 (Shishkina & Wagner 2016), together with (3.1), we obtain ∼Pr 5/4 Ra. And, indeed, DNS (dashed curves in figure 4) supports that a constant Pr 5/4 Ra ≈ 9 × 10 10 determines the onset of the detached plumes regime.
At higher Ra, plumes will detach faster, and for sufficiently large Ra one finds multiple plumes detaching from the thermal BL. This phenomenon was reported in Passaggia et al. (2017) for Ra = 9 × 10 14 , where plumes were visible immediately after entering the convectively unstable region.

Oscillations
A laminar flow in SHC can be thought of as a configuration of two convective flows in subcells meeting in the centre and circulating in opposite directions. While talking about oscillations, we refer to a horizontal movement of these two large structures and analyse an oscillatory movement at the location where the two rolls meet. This location oscillates periodically around the geometric centre of the cell and thus breaks its symmetry.
Following the same strategy as in the previous section, the onset of oscillations can be described in a simplified way as follows. Assume there exists a temperature fluctuation in one of the subcells near the centreline, which, due to buoyancy forces, leads to a local velocity change v of a flow parcel travelling upwards (relative to the base flow) against the viscous forces. As the speed increases, the pressure  drops according to Bernoulli's theorem, which consequently initiates a horizontal pressure gradient between the two convective subcells. This essentially reflects the underlying role of the pressure term in the Navier-Stokes equations, which can transfer energy between modes of different directions (Batchelor 1953). Therefore, a vertically directed force (buoyancy) can induce horizontal oscillations. The characteristic velocity of this Stokes-type flow is v ∼ (αg L 2 )/ν and the time scale is T p ≈ L/( v) = ν/(αg L). The stabilizing antagonist here is the viscous force, which acts as preserving for the symmetric flow profiles and the viscous time scale is T ν = L 2 /ν. The oscillations happen as soon as the shear time scale becomes large, unable to smooth out the asymmetric flow profiles at a certain constant value of T ν /T p : T ν /T p = Ra/Pr = Gr. (3.2) Our DNS results ( figure 4a,b) indicate that the critical value is Ra/Pr ≈ 5 × 10 9 , supporting the above described physical picture. This is consistent with the results of Paparella & Young (2002), who found that the transition to a time-dependent flow occurs at Ra/Pr ≈ 1.6 × 10 8 . The discrepancy in the prefactors is explained by different BCs and that their simulations were 2-D. Two remarks about figure 4 should be made. First, especially for low Pr, periodic oscillations exist only near the onset of the instability (figure 4e). With increasing Ra, the flow becomes chaotic (figure 4f ). Second, we cannot identify the regime of oscillations in CHC, but found the onset to a time-dependent and not plume-determined flow with a similar trend. The different plots in figure 4(c-f ) show the time signals of the vertical heat flux, averaged over the heated plates. As discussed, low-Pr flows show oscillations and chaotic behaviour, while for large Pr we find the presence of plumes and a combination of plumes and oscillations. In general, the frequency of detaching plumes is an order of magnitude larger than the oscillatory frequency (see caption of figure 4). It remains to be noted that the locations of onsets to time-dependent flows shown in figure 4(a,b) coincide with Nu and Re transitions as seen in figure 2(a,b).

Dissipation rates
To study how the transition to a time-dependent flow can affect the global scalings (analysed in § 3.1), we now analyse the kinetic ( u ) and thermal ( θ ) dissipation rates and assess the results in the context of the SGL model. Following Ng et al. (2015), we decompose the dissipation rates into their mean and fluctuating parts: This will give us a qualitative understanding about the role that fluctuations have on the mixing process. Additionally, we consider the volume averages restricted to the domain part above the heated plate, · + , where we expect the most turbulent fluctuations. 892 R1-8 For u V one can expect either the BL scalings ∼Re 2 , ∼Re 5/2 or bulk scaling ∼Re 3 . Figure 5 shows a non-monotonic behaviour in all cases. First, for the lowest Re, the scaling shows approximately u V ∼ Re 2 behaviour, which corresponds to regime I * l with Nu ∼ Ra 1/4 . As Re increases, we observe a rather rapid increase of u V leading to positive slopes in the compensated plot. The sudden increase in u V is accompanied by a region where the dissipation of the mean flow starts to drop. The total kinetic dissipation rate in the region above the heated plate u + increases even stronger and, for high Ra, most of the energy dissipates inside this region. The value of Re where the first dissipation increase occurs correlates strongly with the transition to a time-dependent flow, as indicated by the vertical dashed lines. Subsequently, the curves drop again to slopes in between ∼Re 5/2 and ∼Re 3 . For high Re and especially for low Pr, we observe that the contribution from the mean flow u V is no longer dominant, which matches the observations of Mullarney et al. (2004) and Scotti & White (2011) that turbulent fluctuations start to become dominating in HC. For our highest Ra and Pr = 1 (figure 5b), a transition to a turbulent regime ∼Re 3 appears, but more data points at higher Re are needed to extend this trend. Another observation one can make from figure 5(a,d) is that, for increasing Re, u V and u + first converge and then slightly diverge again. This is explained by the fact that the region where a turbulent flow is present starts above the heated plate, but then spreads over an increasingly larger volume of the domain.
In figure 6 the thermal dissipation rate is analysed in a similar way. Other than for the kinetic dissipation, there is no observable effect from the onset of the instabilities. Moreover, it is evident that the contributions of turbulent fluctuations is small for all studied Ra and that the total thermal dissipation is well described by its mean field contribution. Only for our largest Ra and only above the heated plate does the mean flow dissipation deviate slightly from the total thermal dissipation. The scaling is approximately θ ∼ Re 1/2 to θ ∼ Re 3/4 and nearly constant.
In summary, we found a strong enhancement of u in the vicinity of the onset of the first instabilities. These locally occurring changes can cause the sharp scaling transitions, as observed in § 3.1, suggesting not a 'scale-free' region. The contributions of the mean dissipation u gradually decrease, and for large Re we observe u ∼ Re 3 , which hints towards a transition to a turbulent regime. The temperature fluctuations, for all studied Ra, contribute little to the total thermal dissipation rate, in contrast to the situation of the kinetic dissipation.

Conclusions
Long-runtime DNS were conducted for several decades of Ra and Pr = 0.1, 1 and 10, for classical and symmetrical HC, in order to investigate the global scaling relations and the flow dynamics. The obtained results can be summarized as follows.
First, for the same parameters (Ra, Pr), the SHC and CHC systems provide nearly the same heat and momentum transport (Nu, Re). Thus, we conclude that SHC set-ups can serve as a good alternative to CHC in studying HC systems, which may give a better experimental accuracy, since it gets along without isolating the critical hot wall. The Nu versus Ra scaling analysis for both set-ups showed evidence for regimes I l and I * l , according to , as found previously in Shishkina & Wagner (2016) and Ramme & Hansen (2019). Further, the Nu evolution suggests another transition phase for Ra > 10 10 and Pr 1, which we found to be presumably related to the transition from a steady to a time-dependent bulk flow. For our highest Ra = 3 × 10 12 , both Pr show a slope of Nu ∼ Ra 0.24 .
Second, the analysis of the dynamics of HC systems reveals three different unsteady flow regimes: detached plume regime, oscillatory regime in SHC and chaotic regime. The onsets of the former two instabilities have been obtained theoretically up to a constant and were confirmed by our DNS data. Detaching plumes dominate high-Pr flows and are found above a critical Ra Pr 5/4 ≈ 9 × 10 10 , while the oscillatory instability starts at a Ra/Pr ≈ 5 × 10 9 and is therefore dominating especially in small-Pr fluids. A subsequent examination of the kinetic and thermal dissipation rates showed that the onsets of these instabilities coincide with a strong increase in the total kinetic dissipation and a simultaneous decrease in its mean field contribution. Our DNS show also that velocity fluctuations become the dominating part of u V , while the temperature fluctuations contribute only a little to θ V (less than 5 %).
Further experimental or numerical investigations for Ra > 10 12 are absolutely crucial for verifying of the other regimes of the SGL model and for the understanding of the role of buoyancy forcing on the ocean dynamics.