Weakly nonlinear bifurcation behaviour at the onset of instability in horizontal convection using an advective buoyancy approximation

Horizontal free convection is studied under an extended buoyancy approximation constructed to capture centrifugal effects in strongly rotating regions within buoyancy driven flows. Under this approximation, the Gay-Lussac parameter ( Ga ) arises in the momentum equation to characterise non-Boussinesq behaviour. The scaling of average Nusselt numbers ( Nu avg ) are determined, as is its dependence on Ga . Increasing Ga is found to decrease Nu avg when convection dominates the flow, while no influence is detected in the low Rayleigh number regime dominated by conduction. The influence of Ga on the onset of the time-dependent regime is investigated and the critical Rayleigh number corresponding to different Ga are determined. Results demonstrate that increasing Ga stabilises the flow, delaying transition to the time-dependent regime. An Orr – Sommerfeld type stability analysis is performed to determine the local stability at different horizontal stations. A zone of convective instability is first detected, growing from the hot end of the enclosure, at Rayleigh numbers three orders of magnitude lower than the predicted global stability threshold at the maximum investigated Ga value. Ga is found to play a key role in determining the preferred orientation of instability roll structures, with increasing Ga being accompanied by a bias toward transverse-roll instabilities over longitudinal rolls. The Stuart — Landau model is used to reveal that the non-linear characteristics of this unsteady transition are consistent with a supercritical Hopf bifurcation.


Introduction
The Boussinesq approximation (also known as the Oberbeck-Boussinesq approximation; OB) is almost universally adopted to model buoyancy driven fluid flows.It is predicated on the insight that variations in fluid density may be neglected except within the gravitational (buoyancy) term of the momentum equation.This strictly holds only for small temperature gradients that confine its validity to small ranges of scientific and industrial applications.In recent decades, several remedies have been proposed to improve the accuracy of the OB approximation within its incompressible framework over a larger domain of applicability.A review of the non-OB approximations and formulations adopted within each category may be found in Ref. 1.
The set of non-OB approximations may be divided into two groups: one based on a compressible framework, the other on an incompressible framework.The Gay-Lussac approach is built upon the incompressible equations, extending consideration of density variations to momentum terms beyond the gravitational buoyancy term [2][3][4].The Gay-Lussac parameter Ga = βΔθ emerges from this treatment, characterising the departure from the conventional Boussinesq approximation through a factor (1 − GaΘ) modifying the respective advection and convection terms of the momentum and energy equations.Ref. [5] proposed an approximation in which density variations were extended solely to the momentum advection term.This captures centrifugal effects in rapidly rotating regions that are absent from the conventional Boussinesq approximation [5][6][7].In this study, this advective approximation is applied for local stability analysis of a special type of free convection known as horizontal convection [8].Horizontal convection (HC) describes free convection due to non-uniform buoyancy supply across a horizontal boundary, and its study is motivated by applications to Earth's ocean currents [9], planetary mantles [10] and industrial processes such as glass melting [11].
A limited number of studies have been conducted to study HC under the OB approximation.Hossain et al. [12] explored small aspect ratios of height to length down to H/L = 0.001 up to Ra = 10 16 at Pr = 6.14.Their findings revealed the importance of enclosure depth in dictating the flow dynamics, and provided evidence that end-wall effects swamp the enclosure when H/L > 0.1.Below this, heating and cooling occur only within approximately 4H of each end wall.Mayeli & Sheard [13] International Communications in Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ichmthttps://doi.org/10.1016/j.icheatmasstransfer.2023.107134used an entropy generation measure under a Gay-Lussac type approximation to distinguish flows dominated by either convection or conduction.They reported that when Ga≳0. 5 and Ra ≈ 6 × 10 5 , these processes were balanced.Tsai et al. [14] examined the non-linear evolution of the Hopf instability breaking the time-invariant symmetry of two-dimensional horizontal convection, finding the bifurcation to be supercritical.They considered different horizontal temperature profiles, finding a step-shaped profile to be more unstable compared to smoother variations.Passaggia et al. [15] also considered a step-shaped temperature profile at Pr = 1 in an enclosure with H/L = 0.25.They used a global linear stability analysis to determine critical Rayleigh numbers of 2 × 10 7 and 1.7 × 10 8 for no-slip and free-slip boundary conditions, respectively.
Tsai et al. [16] recognised that laminar two-dimensional HC first destabilizes within its convective horizontal boundary layer.They conducted an Orr-Sommerfeld type stability analysis for a HC system with an aspect ratio of 0.16 and for Prandtl numbers spanning 0.1 ≤ Pr ≤ 10, finding that the flow remains locally stable up to Ra = 5 × 10 7 .Beyond this, a two-dimensional transverse-roll instability emerges, which gives way to a subsequent three-dimensional instability comprising streamwise longitudinal roll vortices consistent with the three-dimensional simulations undertaken in Gayen et al. [17].
Recently, Mayeli et al. [18] conducted a global linear stability analysis of HC at a unit Prandtl number under the advective buoyancy approximation (referred to as centrifugal buoyancy approximation in that study) concluding that across the physical range of Ga, the flows are unstable to infinitesimal oscillatory three-dimensional perturbations beyond Ra ≥ 4.23 × 10 8 .
Presently, our understanding of horizontal convection under the Gay-Lussac framework is incomplete.This study seeks to close several remaining knowledge gaps.Following the methodology described in § 2, verification of the consistency of the centrifugal buoyancy formulation is provided in § 3. § 4 will elucidate the two-dimensional kinematics via vorticity fields and the fields capturing the non-OB effects.§ 5 reveals the influence of Ga on the Nu avg -Ra relationships.§ 6 details a local Orr-Sommerfeld type stability analysis.§ 7 elucidates linear and nonlinear aspects of the two-dimensional Hopf instability breaking the steady solution branch, while § 8 elucidates the corresponding aspects of the first three-dimensional instability.Conclusions are drawn in § 9.

Methodology
Following [6,7], the effect of density variations are considered within the momentum advection term in addition to the gravity term considered under the conventional OB approximation.The continuity, momentum and energy equations may be written as ( As with the OB approximation, here viscous heat dissipation is neglected, thermal diffusivity and kinematic viscosity are assumed constant, and a linear state equation relating density to temperature is assumed, ρ/ρ 0 = 1 − βθ.The equations may be non-dimensionalised by introducing the following scalings, The dimensionless equations then reduce to In Eq. ( 6), Θ is the dimensionless relative temperature.Observe the emergence of the Rayleigh number Ra = gβΔθL ref 3 /να, describing the relationship between buoyancy and dissipation, and the Prandtl number Pr = ν/α, describing the relationship between viscous and thermal dissipation.Eq. ( 6) differs from the conventional OB approximation through the additional term on the right-hand side, which acts on strongly accelerating regions in the fluid.The strength of these deviations is described by Gay-Lussac parameter Ga, with Ga→0 recovering the classical OB approximation.Ga also appears in the linear density state equation, and thus for a physical density ratio ρ/ρ 0 > 0, regarding the dimensionless temperature bounds ( − 0.5 ≤ Θ ≤ 0.5), the allowable range for the Gay-Lussac parameter is 0 ≤ Ga ≤ 2. Numerical results are computed up to Ga = 2 as a ceiling on its allowable range.This work considers horizontal convection in a rectangular enclosure with aspect ratio H/L = 0.16, under the centrifugal buoyancy approximation.The Prandtl number is fixed at Pr = 1, representative of gases such as steam or air.The system is shown in Fig. 1a.All boundaries permit no fluid slippage (U = 0).Side and top boundaries are adiabatic (∂Θ/∂n = 0), while the bottom boundary features a linear temperature profile.This condition creates a buoyancy imbalance along the bottom boundary, compelling the fluid to circulate in an anticlockwise direction for any non-zero Rayleigh number [19].Eqs.

Nomenclature
(5-7) are numerically solved using an in-house code employing a nodal spectral element method for spatial discretisation, and a thirdorder operator-splitting time integration scheme built upon backwards differentiation.The code has been used and validated across numerous natural convection studies [12][13][14]16,18].The mesh discretising the flow domain (Fig. 1b) is constructed with elements concentrated toward both bottom boundary and right side-wall, where the HC boundary layer and side-wall plume structures are expected to be strongest.Mesh independence was established previously for this model in Ref. 18.To summarise, mesh independence was sought at the ceiling of the explored Ra − Ga parameter space, namely Ra = 4 × 10 8 and Ga = 2.Both hrefinement (variation in the number of spectral elements in the mesh) and p-refinement (variation in the polynomial order of the elements) were conducted.It was determined that a mesh composed of 20 spectral elements having a polynomial order of 30 was able to reproduce an integral norm of the velocity field to within 2 × 10 − 7 % of higherresolution cases.This configuration is used throughout the present study.

Verification of the consistency of the centrifugal buoyancy formulation
The formulation employed in this study is specifically the generalised formulation for centrifugal buoyancy in the inertial frame proposed in Lopez et al. [5] (their eq.2.12).This section provides evidence to demonstrate that this formulation is a consistent approximation of the full governing equations, using the formalism developed by Gray & Giorgini [20].
In Gray & Giorgini, the starting point was the governing continuity, momentum and energy equations for a Newtonian fluid with variable properties.Linear approximations capturing the pressure and tempera-ture dependence of density, specific heat capacity, viscosity, thermal expansion and thermal conductivity were employed.Following nondimensionalisation, the governing equations were expressed containing dimensionless constants ε 1 , …, ε 11 .The authors explained that the traditional Boussinesq approximation is recovered by neglecting the terms multiplying these constants on the basis that they are very small.The first of these constants, ε 1 , is the product of the volumetric thermal expansion coefficient and the characteristic temperature difference, or Ga in the present work).In Gray & Giorgini's derivation, terms across the governing equations that are pre-multiplied solely by ε 1 are listed in Table 1, along with their equivalents under the present notation.
The formulation employed herein arises if the centrifugal buoyancy term Θ (U⋅∇)U is sufficiently large that it cannot be neglected, despite being pre-multiplied by ε 1 (Ga).To justify the retention of this term alone, it must be demonstrated that its magnitude is significantly larger than that of the other ε 1 terms.
To verify this, solutions were computed at Ga = 0 and Ra = 4 × 10 5 , 4 × 10 6 , 4 × 10 7 and Ra = 4 × 10 8 .The magnitude of each term in Table 1 was evaluated and the maximum values were recorded over time.The results are plotted in Fig. 2, with the four largest terms being visible.It is evident that the centrifugal buoyancy term Θ(U⋅∇)U is significantly stronger than the other terms (with the exception of a brief time early in the transient startup phase at the lowest Rayleigh number, when the Θ∂U/∂t term is prominent).
Table 2 summarises the maximum magnitudes at steady state for each term in Fig. 2. The centrifugal buoyancy term exceeds the next largest terms by between one and four orders of magnitude across these Rayleigh numbers.This confirms that situations may exist where the centrifugal buoyancy term is sufficiently large that it contributes to the flow dynamics, while all other ε 1 terms remain negligible.Of particular importance, note that the terms arising from the mass conservation equation, ∂Θ/∂t, (U⋅∇)Θ and Θ∇⋅U, are all much smaller than the centrifugal buoyancy term, confirming that their omission does not violate mass conservation under the approximation used herein.

Dependence of flow structure on Ga
The effect of Ga on the temperature field was already investigated in Ref. 18.In this study, its effect on the flow fields across three orders of magnitude beyond Ra = 4 × 10 5 is elucidated.In Fig. 3, solutions are visualized via their vorticity fields.Complementing Fig. 3, the magnitude of the Gay-Lussac modification to advection, |GaΘ((U • ∇)U ) |, is displayed in Fig. 4.
In horizontal convection, the fluid experiences density variations due to imposed temperature differences imposed from the horizontal boundary supplying buoyancy (here the bottom boundary).As a result, regions with different densities develop, and the fluid begins to circulate to equalize these density differences.This process typically leads to the formation of an enclosure-scale overturning circulation captured here by the vorticity fields plotted at different Rayleigh numbers in Fig. 3.

Table 1
Terms pre-multiplied by ε 1 in the working contained within Gray & Giorgini [20], along with their equivalents under the present notation.Under Gray & Giorgini's notation, tildes (~) denote non-dimensional quantities, subscript "0" denotes reference values, and x, t, T and V are the spatial coordinate, time, temperature and velocity vector.

Mayeli and G.J. Sheard
At Ra = 4 × 10 5 (Fig. 3a), there is a main vorticity core that has more strength (and consequently stronger circulation) close to the cooling section over the left half of the enclosure.This is the region where warmer fluid cools and descends to complete the overturning circulation.This core becomes stronger and compressed toward the cooling section as Ga is increased.The difference between the vorticity fields at the top frame of this figure is of order O ( 10 2 ) and has almost a symmetric distribution with larger difference along the four sides except for four corners.The non-OB effects at this Rayleigh number (Fig. 4a) are active in the cooling section and they become stronger with increasing Ga, such that the non-OB field at Ga = 2 is larger by an order of magnitude than at Ga = 0.5 across the bottom-left region of the enclosure.By increasing the Rayleigh number to Ra = 4 × 10 6 (Fig. 3b), there is a similar but stronger overturning circulation that extends over the physical domain.Similar to the previous state, by increasing Ga, this circulation becomes concentrated toward the cooling section with a slightly more strength in its central region.At this Rayleigh number, the difference has a symmetric distribution of order O ( 10 3 ) along the horizontal axis in the middle of the physical domain and the difference primarily arises from contraction of the main core toward the cooling section with increasing Ga.The non-OB effects at Ra = 4 × 10 6 (Fig. 4b) distinct from Ra = 4 × 10 5 are active in both of the cooling and heating sections and they become stronger by increasing Ga.At both Ra = 4 × 10 5 and Ra = 4 × 10 6 , the leftward shift of the main circulating core with increasing Ga is accompanied by the formation of a weak counterclockwise rotation at Ga = 1.5 and 2 over the heating section, which is completely absent at Ga ≤ 1.
At Ra = 4 × 10 7 (Fig. 3c), there is an overturning clockwise circulation in the plume region at the bottom-right region of the enclosure at Ga = 0.5, 1 and 1.5 and a strong counter-clockwise circulation over the bottom side due to thermal boundary layer.At this Rayleigh number, the ascending plume at the right is strong so that the interaction of the fluid circulation with the horizontal and vertical sides creates two circulations across the top-right corner of the enclosure.Excitingly, the non-OB effects at this Rayleigh number presented in Fig. 4c indicate that, by increasing Ga up to unity, the non-OB effects are augmented, but beyond that, these effects recede.The net non-OB effects at Ga = 2 causes shifting of the main core to the left with a lower strength.Similar to two previous cases, a weak counter-clockwise circulation is formed over the heating section whose power is so weak compared to the main one.The difference of the vorticity fields at Ga = 0 and Ga = 2 at this Rayleigh number is of order O ( ) and as it can be seen, it mainly comes from shifting the main core leftward at Ga = 2 compared to Ga = 0.The non-OB effects at this Rayleigh number do not follow the previous analysis.
Finally, at Ra = 4 × 10 8 , an overturning plume in the heating region is clearly visible.The flow pattern is similar to Ra = 4 × 10 7 but here the fluid has sufficient buoyancy to reach to the top of the enclosure and form an overturning circulation at the top-right zone.As seen, the strength of the circulations over the plume region at both Ra = 4 × 10 7 and 4 × 10 8 has weakened by increasing Ga.At Ga = 0.5, 1 and 1.5, the plume interaction with the vertical and horizontal walls has created the two counter-clockwise vortices attached to the surfaces across the topright region but at Ga = 2, these vortices are annihilated.The difference of the vorticity fields at this Rayleigh number is of order O ( 10 5 ) and it mainly happens across X≳0.75.The non-OB effects at this Rayleigh number increases with increasing Ga up to Ga = 1.5 but weaker non-OB effects are observed thereafter.

Dependence of Nu avg on Ga
It is apparent from the previous section that these flows can be visibly altered by varying Ga, particularly across Rayleigh numbers producing convective behaviour.The corresponding effects on heat transfer are now considered.The average Nusselt number is computed at Ga = 0, 0.5, 1, 1.5 and 2 for Rayleigh numbers up to Ra = 5 × 10 8 .Results are shown in Fig. 5. Nu avg is seen to remain almost constant up to Ra ≈ 10 5 .This is a known behaviour [21], but it is now shown to persist for nonzero Ga.Beyond Ra ≈ 10 6 , Nu avg increases as the flow enters the convection-dominated regime, and a dependence on Ga also emerges as the data sets depart from each other, with larger Ga depressing Nu avg .The deviations develop at Rayleigh numbers approximately an order of magnitude beyond the Rayleigh number where the Nusselt numbers first begin to rise from the low-Ra baseline.Calculations show around 17% difference in Nu avg between Ga = 0 and 2 at Ra = 5 × 10 8 .The trends are highly linear on the log-log axes, indicating a power-law dependence.At Ga = 0, the average Nusselt number beyond Ra ≈ 5 × 10 6 is found to adhere to Nu avg ∼ Ra 0.206 , while for Ga = 2 presents Nu avg ∼ Ra 0.257 .These exponents are respectively close to 1/5 and 1/4, which are notable power law scalings theorised for HC flow [22,23].
The power-law exponents for the Nu avg -Ra relationships in the convective regime for each Ga are given in Table 3.The scaling at Ga = 0 is consistent with the Rossby's 1/5 scaling for the average Nusselt number [24] and the exponent monotonically rise to ∼ 1/4 with increasing Ga.Lower Nu avg values at higher Ga may be attributed to the local Nusselt number behaviour [18].In the conduction-dominated regime, Ra = 10 5 for example (Fig. 3b), Nu loc along the bottom surface is symmetric and there is almost no difference among different Ga.However, in the convection-dominated regime, Ra = 4 × 10 8 for instance (Fig. 3c), Nu loc is decreased by increasing Ga value except for a small range of 0.45≲X≲0.75.
A possible explanation for the reduction in Nusselt number with increasing Ga may be developed from Eq. ( 6), in which the advection term on the LHS is eroded by positive Θ at non-zero Ga due to the last term on the RHS (i.e.effectively there is a 1 − GaΘ prefactor on the advection term).As velocities and gradients are strongest in the boundary layer at the hotter end of a horizontal convection flow, increasing Ga may act to dampen the action of advection in this region, which in turn would reduce thermal gradients, reducing the Nusselt number.

Local stability under the centrifugal buoyancy approximation
In a recent global linear stability analysis of this system [18], the predicted instability modes were predominantly located in the vicinity of the ascending plume at the hot end of the enclosure.However, the literature provides evidence that instability may be seeded by the convective boundary layer adjacent to the horizontal boundary supplying buoyancy to the flow [17].To address the question as to whether the horizontal convection boundary layer is the source of instability, and  the role, if any, of non-OB effects, consideration is now given to the local stability of the two-dimensional base flows.The basis of this analysis is to extract one-dimensional flow profiles at discrete horizontal stations along the breadth of the enclosure.A key underpinning assumption of this analysis is that the flow is approximately horizontal (i.e.∂/∂x, ∂/∂z≪∂/∂y).This permits the problem to be constructed as an Orr-Sommerfeld type stability analysis to reveal the local stability at each station to small three-dimensional time-dependent perturbations.To derive the eigenvalue problem, a base flow containing perturbations is defined in the general form φ = φ b (y) + εφ(y)e i(αX+βZ− ωt) , (9) in which φ is any of the flow variables, ε is a small constant and φ is a complex eigenfunction.The perturbation comprises travelling wavenumbers α and β in the respective X and Z directions, and a frequency and growth rate described by ω.Here, the flow under investigation is not strictly parallel, thus, we restrict local stability analysis to the regions where the one-dimensional flow assumption is approximately valid.This region is found by comparing the absolute ratio of the base flow velocity components.In general, regions where the ratio of vertical to horizontal velocity is below 10% are considered as appropriate regions for this analysis.The same criterion is considered for horizontal velocity gradients as further support for a one-dimensional flow assumption in the region of interest.
The eigenvalue problem is derived by substituting Eq. ( 9) into the governing eqs.(5-7), and retaining terms of order O(ε).The resulting equations are where D and ( ′ ) represent partial derivatives and differentiation with respect to the vertical direction, respectively.The wave number is expressed by k 2 = α 2 + β 2 .Eqs. (10− 11) yield an eigenvalue problem in a general form of Ax k = ω k Bx k , where eigenvector x k contains the perturbation Ṽ and θ coefficients, while matrices A and B contain the respective left and right side of Eqs.(10)(11).The complex eigenvalue ω = ω Re + iσ Im evolves as exp(σ Im )[cos(ω Re ) + isin(ω Re ) ]. Hence the imaginary part σ Im expresses the exponential growth rate of the instability, and the real part ω Re yields its angular frequency.Typically, perturbations are considered that are invariant in either the streamwise (X) direction or the spanwise (Z) direction [16].These are conventionally termed longitudinal-roll or transverse-roll instabilities, respectively, due to the orientation of perturbation vortices they describe.They are respectively investigated by setting α = 0 or β = 0 and seeking the value of the alternate wavenumber maximising σ Im .
A pseudo-spectral method is employed to solve this eigenvalue problem, following Ref.[16] and references therein.Accurate performance of the solver was verified by computing the well-established critical eigenmode in plane Poiseuille flow [25].The exact critical Reynolds number Re cr = 5772.22and wavenumber α cr = 1.02056 were obtained.An additional test was conducted on a buoyancy-driven problem, Rayleigh-Bénard convection, where the expected critical Rayleigh number Ra cr = 1707.76and wavenumber α cr = 3.177 [26] were obtained.
To initiate the stability analysis, it is necessary to first determine the region where the one-dimensionality assumption is valid.The parallel flow assumption is tested using measures that reduce to zero as parallel flow is achieved for the extracted base flow velocity profiles at Ra = 10 9 in Fig. 6a and b  Regarding the eigenvalue problem, a mesh independency test is performed for the number of quadrature points required and it is found that Chebyshev polynomial basis function of order 30 are sufficient to keep the stability results independent of further resolution.The horizontal velocity and temperature profiles at an arbitrary location from the considered range (here X = 0.8 for example) are plotted in Fig. 6c  and d.These profiles serve as the base flow shown with an overbar in Eqs. 10 and 11.
The first step is to determine where and when the flow first becomes locally unstable.In this Fig. 7, growth rates corresponding to the transverse (β = 0) and longitudinal (α = 0) roll instabilities are plotted from stations X = 0.70, 0.75, 80 and 0.85, and Ga = 0, 1 and 2. Local convective instability is indicated by positive growth rates.These results show that the flow becomes more unstable as the horizontal station is moved to the right, and increasing both Ga and Ra also promotes an elevation in growth rate.Up to Ra = 5 × 10 6 , the flow is locally stable to both transverse and longitudinal roll instabilities at Ga = 0 and 1 (Fig. 7a and b) over the considered length but the growth rate of the perturbations at Ga = 2 and X = 0.85 (Fig. 7c) exceeds σ = 0 at wavenumber k = 19 for the first time.Linear stability analysis for this problem [18] indicates that flow becomes globally unstable at Ra = 4.23 × 10 8 and Ga = 2. Fig. 7c indicates that the flow becomes locally unstable two orders of magnitude lower than the reported critical Rayleigh number for the global stability analysis at Ga = 2.At this Rayleigh number, there is no preference between the transverse and longitudinal roll instabilities in crossing σ = 0 but there is clear difference between the two roll types at higher Rayleigh numbers.
Figs. 7d-f present local stability analysis results at one order of magnitude higher Rayleigh number, i.e.Ra = 5 × 10 7 .As seen, at this (and higher) Rayleigh numbers, perturbations growth rates corresponding to the transverse rolls are higher than the longitudinal rolls in the σ − k space, which indicates the precedence of the transverse rolls instabilities compared to the longitudinal ones.The data in Fig. 7d also indicates that at Ra = 5 × 10 7 and Ga = 0, transverse roll instabilities become neutrally stable at X = 0.85.For comparison, the critical Rayleigh number at Ga = 0 predicted by the global linear stability analysis of Ref. 18 was an order of magnitude higher at Ra cr = 6.46 × 10 8 .Presented local stability results at the same Rayleigh number at Ga = 1 (Fig. 7e) reveals that a band of longitudinal roll instability wavenumbers have become unstable (σ > 0), while a wider band of wavenumbers corresponding to the transverse roll instability have already passed this threshold.Similar to Ga = 0, this shows one order of magnitude lower value for the critical Rayleigh number at Ga = 1 compared to the predicted critical Rayleigh number from the global linear stability analysis [18].
As previously mentioned, global linear stability analysis for this problem under the centrifugal approximation predicted critical Rayleigh numbers over Ra ≥ 4.23 × 10 8 for three-dimensional perturbations over all Ga.Fig. 8 shows the local stability analysis results conducted at selected Rayleigh numbers and their corresponding critical Ga predicted by global linear stability analysis [18].The Rayleigh numbers are Ra = 4.23 × 10 8 ,5.5 × 10 8 and 6.46 × 10 8 that are neutrally stable at Ga = 2, 1 and 0, respectively.As seen, at Ga = 2 the maximum of the perturbation growth rates at all X-locations exceed the neutral stability threshold in Figs.7c,f,i, but the optimum of the perturbation growth rates at X = 0.7 remains in the stable region (σ < 0) for Ga = 0 and 1.A feature of the data in Fig. 8 is the location of the maximum wavenumber of the transverse and longitudinal rolls at high Rayleigh numbers.Results indicate that by increasing Ga in the convection-dominated regime, transverse roll are preceded to unstable region but the maximum of the perturbation growth rates occur at a higher wavenumber compared to the longitudinal rolls.Local stability analysis results also show a larger dominant wavenumber for the transverse rolls with increasing Rayleigh number.A reduction in distance between transverse rolls was observed in three-dimensional simulations of HC performed by Gayen et al. [17], consistent with the present analysis.
Marginal stability was determined across the parameter space, and the resulting curves are plotted in Fig. 9 for different X-locations and Ga = 0, 1 and 2. Transverse rolls are found to be slightly more unstable than longitudinal rolls.At Ga = 0 and 1, precedence of the transverse rolls is clear against the longitudinal ones at X = 0.85 and lower X-locations, but at Ga = 2, transition of the transverse and longitudinal instability rolls to the unstable zone occurs together and remains unaltered up to Ra≲9 × 10 7 .It is also found that, for Ra≳10 9 (especilly for Ga = 1 and 2), the transverse rolls from a relatively more stable region farther from the hot end-wall corner exceed the longitudinal instability rolls from a location closer to the right side wall, which reflects the strength of the non-Boussinesq effects in the fully convection-dominated regime.This is in agreement with [17], where transverse roll features in the HC boundary layer were overwhelmed by longitudinal rolls.

2D instability: Insights from the nonlinear Stuart-Landau model
In this section, nonlinear simulations are used to investigate the transition to time-dependent two-dimensional flow.A useful tool is the Stuart-Landau equation [27,28] which describes the growth and saturation of a complex amplitude measure of an instability A(t), where σ and ω are respectively the amplitude growth rate and angular frequency at asymptotically small amplitudes.Coefficients l and c serve to describe the weakly nonlinear departure of the mode evolution from the linear regime.Higher-order terms are truncated from the right hand side of Eq. ( 9).If A(t) takes the form |A|exp(iϕ) with magnitude |A| and phase angle ϕ, it is possible to decompose Eq. ( 9) into real and imaginary parts Therefore, if we record A(t) for a growing or decaying instability mode, insight may be gleaned by considering Eq. ( 13), specifically the relationship between d(log|A|)/dt and |A| 2 .At |A| 2 = 0, d(log|A|)/dt expresses the linear growth rate σ, while its gradient yields − l.The sign of this constant reveals the weakly non-linear characteristics of a developing instability mode, where positive and negative slopes (− l > 0 and − l < 0) respectively describe subcritical and supercritical bifurcations.A supercritical bifurcation permits only a stable branch below the critical parameter, whereas a subcritical bifurcation permits the unstable branch below the critical parameter.This technique has been extensively used for external flows over bluff-bodies such as toroidal bodies [29], cylinder with square cross-section flow [30] and wakes behind a sphere [31].
Here the envelope of the time-periodic fluctuation of the average Nusselt number is used as the amplitude measure.The procedure used to determine the stability threshold is as follows: first a time-periodic solution at an unstable Rayleigh number is obtained for each considered Ga value.Each of these solutions forms the initial conditions for several subsequent simulations conducted at several incrementally smaller Rayleigh numbers.Nusselt number time histories are obtained until the flows stabilise.The envelopes are then determined, giving |A(t) |, from which the growth rates may be found.Fig. 10 shows a representative case of Stuart-Landau analysis at Ga = 0.5.In this case, Ra = 3 × 10 9 yields a time-periodic state which is used as an initial condition.An evolution of the average Nusselt at Ra = 2.75 × 10 9 is plotted in Fig. 10a.Fig. 10b shows the corresponding plot of d(log|A|)/dt against |A| 2 that reveals the decay rate of σ as |A| 2 →0.The process is repeated for other Rayleigh numbers and the decay rate is calculated for each case and plotted in Fig. 10c.A linear extrapolation of the obtained data to zero growth rate provides an estimation of the critical Rayleigh number equal to Ra cr = 2.87 × 10 9 at Ga = 0.5.
The critical Rayleigh number for other Ga are calculated similarly, and the results are reported in Table 4.As seen, by increasing Ga, the critical Rayleigh number is also increased.In other words, it is concluded that the non-OB effects captured by Ga act to stabilise the buoyancy driven flow transition to the time-periodic regime.
Another means of estimating the critical Rayleigh number is to record the amplitude of oscillation in the Nusselt number time histories after the instability has saturated at several Rayleigh numbers beyond the threshold, and extrapolating these back to zero amplitude.This was conducted, and the results are compared in Table 4.The close agreement between these methods suggests that finite-amplitude unstable states cannot persist below the critical Rayleigh number.This is consistent with a supercritical bifurcation [27].
Returning to Fig. 10(b), the gradient exhibited in the vicinity of the vertical axis (|A| ≈ 0) is negative.This corresponds to a positive coefficient l > 0, and therefore indicates that this two-dimensional instability is supercritical.This is in agreement with the findings of Tsai et al. [14],    as their study showed the OB unsteady transition in HC to also occur via a supercritical bifurcation.
In the next section, the Stuart-Landau model will be used to probe the nonlinear properties of the three-dimensional instabilities arising in these flows.

3D instability: Insights from the nonlinear Stuart-Landau model
In the recent linear stability analysis conducted for threedimensional disturbances in HC under the centrifugal buoyancy approximation [18], critical Rayleigh numbers were determined beyond Ra ≥ 4.23 × 10 8 across the available range of Ga.The instability consistently led to an oscillatory three-dimensional state, though the non-linear properties of this instability are yet to be determined.In this section, the instabilities predicted in Ref. 18 are subjected to the Stuart-Landau model to elucidate their weakly non-linear properties at Ga = 0, 1 and 2.
The envelope of the domain integral of w-velocity oscillations obtained from direct numerical simulation (DNS) are recorded.The simulations employ a Fourier spectral treatment in the spanwise (Z) direction, and the same 2D nodal spectral element discretisation in the X-Y plane as used earlier in this study.The 3D solutions are then naturally periodic over same prescribed Z-wavelength, and resolution is controlled by specifying the number of Fourier modes employed.The number of Fourier modes required for independency had been determined in Ref. 18 to be 8, which is adopted herein.
Figs. 11a,c,e represent the integral of the absolute spanwise (Z direction) velocity oscillations of three unstable modes including Ra = 5 × 10 8 with Ga = 2, Ra = 6.5 × 10 8 with Ga = 1 and Ra = 7.25 × 10 8 with Ga = 0, respectively.Minimum and maximum of the oscillations are specified by two dashed lines in the integral of absolute w-velocity plots.In all cases, the spanwise wavenumber (stated with k in each figure) is set correspond to the maximum perturbation growth rate obtained by linear stability analysis [18].These evolutions are computed between 12% and 18% beyond the critical Rayleigh numbers at each considered Ga value.
Results from the Stuart-Landau analysis are presented in Figs.11b,d, f.From these figures, extrapolation to |A| 2 = 0 yields growth rate σ.The nearly linear trends having a negative slope toward small |A| 2 is evidence that the instability bifurcation is consistently supercritical.Growth rates are also estimated by separately monitoring the time history of kinetic energy in the Fourier mode associated with the underlying linear instability mode.These growth rates are reported for comparison in Table 5, with close agreement being observed.

Conclusions
Horizontal convection is studied under the advective buoyancy approximation within a rectangular enclosure having a ratio of height to length of 0.16 at a Prandtl number of unity.The average Nusselt number was scaled against the Rayleigh number at different Gay-Lussac parameters, concluding that by increasing Ga, the average Nusselt number decreases in the convection-dominated regime (Ra≳2 × 10 6 ), while no difference is seen at lower Rayleigh number.A one-dimensional stability analysis was conducted under the centrifugal buoyancy approximation to investigate local stability to both the transverse and longitudinal roll instability modes.This analysis indicates that the buoyancy-driven flow becomes locally unstable two orders of magnitude lower compared to global stability analysis at Ga = 2.This difference reduces to one order of magnitude for the lower values Ga = 1 and 0. This may serve to explain why the global linear stability threshold predicted in [18] advanced to a lower Rayleigh number as Ga increased.The local stability analysis results also show that as Rayleigh number is increased, the flow becomes locally unstable to a transverse-roll mode ahead of a longitudinal-roll mode.The effect of Ga on the flow transition to the time-dependent regime was investigated via a truncated Stuart-Landau analysis.Critical Rayleigh number of the horizontal convection system were calculated and reported for different Ga and verified against the third order accuracy in time solver.Results indicate that the flow becomes unsteady for the first time at Ra = 1.95 × 10 9 with Ga = 0.It was found that Ga acts as a stabiliser and delays the transition to the timeperiodic regime.The Stuart-Landau results indicate that the buoyancy-driven flow transition to the time dependent regime in horizontal convection occurs through a supercritical bifurcation.The Stuart-Landau model was also applied to the evolution of global 3D linear instabilities, indicating that the flow transitions occur via a supercritical Hopf bifurcation.

Declaration of Competing Interest
None.

Table 5
Comparison of the perturbation growth rates calculated from DNS and Stuart-Landau (SL) model analysis at three Ga values.

Fig. 1 .
Fig. 1.(a) Schematic diagram of the system under consideration, including dimensions, coordinate system and applied boundary conditions.(b) The mesh employed in this study, having 5 × 4 spectral elements (black lines).Grids (blue lines) are drawn through interior quadrature points in this visualisation.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2 . 2
Fig. 2. Time histories of the leading terms pre-multiplied by ε 1 (Ga) at the Rayleigh numbers as indicated and with Ga = 0.

2 Fig. 3 .
Fig. 3. Comparison of vorticity fields at different Ra and Ga values, as stated.The top frame of each set shows the absolute difference in vorticity between the Ga = 2 and Ga = 0 (OB) cases.Subsequent frames show vorticity fields for increasing Ga.Vorticity is scaled by Ra 1/2 for consistent levels across the figure, and streamlines are overlaid to aid comparison.

Fig. 4 .
Fig. 4. Contour plots showing fields of the magnitude of the non-OB term (|GaΘ(U • ∇)U |) scaled by Rayleigh number, plotted at Ra and Ga values as indicated.

Fig. 5 .
Fig. 5. Average Nusselt number against the Rayleigh number at different Ga.
. Results indicate that in the region 0.4≲X≲0.85both ∑ |U|/ ∑ |V| and ∑ |∂U/∂X|/ ∑ |∂U/∂Y| are below 10%, so the local stability analysis is restricted to this region in this study.

Fig. 6 .
Fig. 6.Verification of the parallel flow assumption at Ra = 10 9 at different Ga.(a) A comparison of vertical and horizontal velocity magnitude summing over the depth of the enclosure.(b) Comparison of the horizontal velocity gradients in the X and Y-directions summing over the depth of the enclosure.(c) Horizontal velocity profile at X = 0.8.(d) Temperature profiles at X = 0.8.The horizontal dashed lines in (a) and (b) represents the 10% level appropriate for 1D stability analysis.

Fig. 10 .
Fig. 10.(a) Time history showing the decay of an oscillation in the average Nusselt number obtained at Ra = 3 × 10 9 when evolved at the lower Rayleigh number Ra = 2.75 × 10 9 .Here Ga = 0.5.(b) A plot of d(log|A|)/dt against |A| 2 (c) Instability decay rate against Ra at Ga = 0.5.The linear fit reveals the critical Rayleigh number where it intersects σ = 0.

Fig. 11 .
Fig. 11.Time histories of ∫ |w|dΩ and time derivative of the amplitude logarithm against the square of the amplitude oscillations for (a, b) Ga = 2 and Ra = 5 × 10 8 with k = 45.51 (c, d) Ga = 1 and Ra = 6.5 × 10 8 with k = 48.25 and (e, f) Ga = 0 and Ra = 7.25 × 10 8 with k = 61.63,demonstrating a supercritical behaviour.The solid black circles in (b, d, f) represent the growth rate predicted by the Stuart-Landau model.

Table 3
Exponents (n) for the power-law scalings Nu avg ∼ Ra n captured in the convection-dominated regime at different Ga.

Table 4
Comparison of critical Rayleigh numbers estimated from the growth rates from the Stuart-Landau modelling and from an extrapolation of the Nusselt number oscillation amplitudes computed beyond the critical Rayleigh number.