2 Formulation of the problem

Consider a horizontally oriented plane Couette flow configuration with the downward force of gravity $g$ in the Cartesian coordinates $(x_{\ast },y_{\ast },z_{\ast })$ . Between the two plates placed at $y_{\ast }=\pm h$ , a Newtonian fluid is sheared and heated from above. We assume that the velocity $\boldsymbol{u}_{\ast }=(u_{\ast },v_{\ast },w_{\ast })$ , temperature $\unicode[STIX]{x1D703}_{\ast }$ and density $\unicode[STIX]{x1D70C}_{\ast }$ satisfy the boundary conditions

and that the fluid has kinematic viscosity $\unicode[STIX]{x1D708}$ , thermal diffusivity $\unicode[STIX]{x1D705}$ , and coefficient of thermal expansion  $\unicode[STIX]{x1D6FE}$ . Here $u_{0},\unicode[STIX]{x1D703}_{0},\unicode[STIX]{x1D70C}_{0},\unicode[STIX]{x1D703}_{1},\unicode[STIX]{x1D70C}_{1}$ are constants. Using the Boussinesq approximation, assuming $\unicode[STIX]{x1D70C}_{0}\gg \unicode[STIX]{x1D70C}_{1}$ , the fluid motion is governed by the Navier–Stokes equations

where $\unicode[STIX]{x1D735}=(\unicode[STIX]{x2202}_{x_{\ast }},\unicode[STIX]{x2202}_{y_{\ast }},\unicode[STIX]{x2202}_{z_{\ast }})$ .

If we scale the variables as $t_{\ast }=(h/u_{0})t$ , $(x_{\ast },y_{\ast },z_{\ast })=h(x,y,z)$ , $p_{\ast }-\unicode[STIX]{x1D70C}_{0}gy_{\ast }=\unicode[STIX]{x1D70C}_{0}u_{0}^{2}p$ , $(u_{\ast },v_{\ast },w_{\ast })=u_{0}(u,v,w)$ , $(\unicode[STIX]{x1D703}_{\ast }-\unicode[STIX]{x1D703}_{0})=\unicode[STIX]{x1D703}_{1}\unicode[STIX]{x1D703}$ , and eliminate the density from the equations, we have the non-dimensional equations

subject to the boundary conditions

Note that $\unicode[STIX]{x1D735}$ now represents $(\unicode[STIX]{x2202}_{x},\unicode[STIX]{x2202}_{y},\unicode[STIX]{x2202}_{z})$ . The system (2.3) is controlled by the Reynolds number, Rayleigh number and Prandtl number, defined by

respectively. The definition of $Ra$ is the same as that used in Clever & Busse (Reference Clever and Busse1992), and thus $Ra<0$ represents the stable stratification. Since we have assumed that the temperature is proportional to the density, the system is also equivalent to that treated in Eaves & Caulfield (Reference Eaves and Caulfield2015). The bulk Richardson number $Ri$ in that paper can be related to the Rayleigh number by

where Eaves & Caulfield (Reference Eaves and Caulfield2015) used $(Re,Pr)=(1000,1)$ .

The linear basic flow profile of (2.3) can be found analytically as $[u,v,w,p,\unicode[STIX]{x1D703}]=[y,0,0,p_{b}(y),y]$ , where $p_{b}(y)=P_{0}+Riy^{2}/2$ , whilst we must of course use a numerical approach to find nonlinear equilibrium solutions. Here we seek non-trivial steady states. In the periodic domain $x\in [0,2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}]$ , $z\in [0,2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FD}]$ , the flow is discretised by truncated Fourier expansions in $x,z$ , and a Chebyshev expansion in  $y$ . After Galerkin projection, we have algebraic equations that can be solved by Newton’s method. The basic numerical algorithm is the same as Nagata (Reference Nagata1990), Clever & Busse (Reference Clever and Busse1992) and a further description of the present code can be found in Deguchi, Hall & Walton (Reference Deguchi, Hall and Walton2013).

In the following sections we shall investigate the asymptotic development of the steady solutions at large  $Re$ . In the subcritical stratified shear flow study, the range of the parameters $(Re,Ri)$ at which transition to turbulence could occur, bounded by the critical parameter $Ri=Ri_{c}(Re)$ say, has been a matter of great concern. The asymptotic analysis of the steady solutions may give, if they can be regarded as a skeleton of turbulence, an important clue for the large Reynolds number development of  $Ri_{c}$ .

Throughout the paper we assume that the solutions have so-called shift–reflect and shift–rotation symmetries

These symmetries were found in the steady wavy instability of the roll convection cell (Clever, Busse & Kelly Reference Clever, Busse and Kelly1977), and when the thermal effect is absent, the symmetries are identical to those discussed in Nagata (Reference Nagata1990).

3 Large-scale structures for small Richardson numbers

Our aim in this section is to confirm the scaling $Ri\sim O(Re^{-2})$ noticed in Eaves & Caulfield (Reference Eaves and Caulfield2015), using the computation of nonlinear solutions and formal large Reynolds number asymptotic analysis. In view of (2.5), it is sufficient to show $Ra\sim O(1)$ ; note that $O(1)$ means $O(Re^{0})$ throughout the present paper. Based on the VWI theory, we shall consider the derivation of the scaling using a formal large Reynolds number asymptotic analysis of the fully three-dimensional Navier–Stokes equations. However, since only the streamwise-independent component of the flow, called the ‘vortex’ in the VWI theory, plays an essential role in deriving the scaling of  $Ra$ , we first consider flows independent of  $x$ (later we shall see that the balance of the buoyancy term in the $x$ -dependent part is less crucial). The roll-cell found in the unstable Rayleigh–Bénard problem belongs to this class, and it is well known (e.g. from Chandrasekhar Reference Chandrasekhar1961) that if we write

with $[\overline{\boldsymbol{u}},\overline{p},\overline{\unicode[STIX]{x1D703}}](y,z)\in \mathbb{R}$ in (2.3), then we obtain $Re$ -independent equations

where $(\overline{u},\overline{v},\overline{w},\overline{\unicode[STIX]{x1D703}})=(\pm 1,0,0,0)$ at $y=\pm 1$ . From (3.2a ) it is now clear that $Ra$ must be an $O(1)$ quantity, to maintain the balance between the viscous and the buoyancy terms. This means that the vortex component in VWI feels the effect of stratification when $Ra\sim O(1)$ and the asymptotic structure is destroyed for larger  $Ra$ .

Figure 1. Bifurcation diagram of stratified plane Couette flow for $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},Pr)=(0.460,1.905,1)$ . Drag $D$ is shown versus Rayleigh number  $Ra$ . The thick dashed curve is the two-dimensional roll-cell solution branch for $Re=5000$ ; the open circles are the same result but for $Re=10\,000$ . The thick solid curve is the three-dimensional solution branch for $Re=5000$ ; the filled circles are the same result but for $Re=10\,000$ . The thin blue dotted curve is the mirror-symmetric solution branch for $Re=5000$ . The bifurcation points for $Re=5000$ are shown by the triangles.

If the fluid is unstably stratified ( $Ra>0$ ), the vortex equations (3.2) can have a non-trivial solution in their own right. The dashed curve in figure 1 is such a two-dimensional solution branch, computed for $(Re,Pr,\unicode[STIX]{x1D6FD})=(5000,1,1.905)$ . Here, the deviation of the drag on the wall $D$ (i.e. $x$ – $z$ average of $u_{y}|_{y=1}$ ) from its laminar value $D=1$ indicates the presence of nonlinearity. Waleffe (Reference Waleffe1997) referred to $\overline{u}$ and $(\overline{v},\overline{w})$ as the streak and roll components, respectively, and here we also adopt that terminology. The roll and temperature components of the two-dimensional solution are exactly the same as the roll-cell solution found in the Rayleigh–Bénard problem bifurcated at a certain critical Rayleigh number; see Chandrasekhar (Reference Chandrasekhar1961). The streak component is passively driven by the roll cell, and its magnitude is measured by the drag  $D$ .

On the other hand, if the fluid is stably stratified ( $Ra<0$ ), the thermal effect typically subtracts energy from the roll component, and thus three-dimensionality is essential for sustaining the vortex component via the Reynolds stress. For the isothermal case, Hall & Smith (Reference Hall and Smith1991) superimposed a monochromatic wave to support the vortex component. In the same spirit, we add the wave part of the flow

to $(u,v,w,p,\unicode[STIX]{x1D703})$ . Here, $\unicode[STIX]{x1D6E5}\ll 1$ describes the amplitude of the wave. More formally, we assume here that the leading-order part in the large $Re$ asymptotic expansion is the sum of (3.1) and (3.3). Substituting the leading-order flow into (2.3) and neglecting the small terms, we have the vortex equations (3.2) and the wave equations

If $Ra\sim O(1)$ , the wave equations (3.4) become the single equation

where $\widetilde{p}_{y}=0$ must be satisfied on the walls.

If $Pr\sim O(1)$ , then from (2.5), we have $Ri\ll 1$ . Therefore, (3.5) reduces to the Rayleigh equation for the streak derived in Hall & Smith (Reference Hall and Smith1991). This equation is inviscid and thus has a singularity at the critical position where $\overline{u}$ vanishes. In the critical layer of thickness $\unicode[STIX]{x1D6FF}=Re^{-1/3}$ around that position, the viscous effect in the wave equation must be retained, and thus we need a separate analysis there. The analysis of the singularity is essential to understand the mechanism by which the small wave drives the roll. If the critical layer is horizontal, from the Frobenius expansion near the critical position we can see that $\widetilde{p}\sim O(\unicode[STIX]{x1D6E5})$ . Then, since $\widetilde{w}\sim O(\unicode[STIX]{x1D6FF}^{-1}\unicode[STIX]{x1D6E5})$ from the $z$ component of (3.4a ), we can see that the Reynolds stress driving the roll, which is $O(\widetilde{w}^{2})$ , is amplified in the critical layer. This stress should balance the viscous term of $Re^{-2}\unicode[STIX]{x1D6FF}^{-1}$ , which is the first derivative of the roll shear stress. (Note that the value of the roll component is continuous across the critical layer and hence what is influenced by the wave is the $O(\unicode[STIX]{x1D6FF})$ smaller roll shear stress.) Thus if $\unicode[STIX]{x1D6E5}=\unicode[STIX]{x1D6FF}^{1/2}Re^{-1}$ , the concentrated Reynolds stress produces a finite stress jump in the roll across the critical position. It is possible to deduce an analytic form of the jump in terms of $\widetilde{p}$ , but since it is identical to that of the VWI theory of Hall & Smith (Reference Hall and Smith1991), we do not repeat any further analysis here. If the Richardson number becomes $O(1)$ , (3.5) has the form of the Taylor–Goldstein equation for the streak; however for a three-dimensional flow, this is only possible when $Pr\sim O(Re^{-2})$ if $Ra\sim O(1)$ . The reason why the wave part less feels the effect of stratification than the vortex part is that the effect of shear advection acting on the wave is much stronger than that acting on the vortex; recall that the vortex component is advected by a small roll velocity component.

A similar tripartite interaction between the roll, streak and wave, has been treated in the SSP theory by Waleffe (Reference Waleffe1997), but the excitation of the wave was explained by the Orr–Sommerfeld equation rather than the Rayleigh equation. Due to the viscosity remaining in the wave equation, the Reynolds number scaling of the flow field is not fully determined in the SSP framework. In this sense, the reduced description of coherent structures proposed by SSP is considered to be an intermediate reduction of the Navier–Stokes equations to the fully asymptotic VWI approach. However one of the advantages of retaining viscosity is that we can avoid any complication associated with the critical layer singularity: see Hall & Horseman (Reference Hall and Horseman1991). Recently, Blackburn et al. (Reference Blackburn, Hall and Sherwin2013) derived a ‘hybrid’ model of coherent structures in plane Couette flow using such a regularisation. Subsequently Beaume et al. (Reference Beaume, Chini, Julien and Knobloch2015) derived a similar model to analyse a forced flow.

If a stably stratified fluid is subjected to a shear, there is a competition between the stabilisation effect due to the stratification and the destabilisation effect due to the wave instability. As shown in Clever & Busse (Reference Clever and Busse1992), this competition leads to a subcritical bifurcation of the three-dimensional solution from the two-dimensional roll-cell solution, and the bifurcated branch can be continued to $Ra<0$ . The solid curve in figure 1 shows the three-dimensional solution branch for $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=(0.460,1.905)$ , corresponding to the ‘narrow geometry’ used in Eaves & Caulfield (Reference Eaves and Caulfield2015). Decreasing $Ra$ from zero, at the other end of the solution branch the flow field acquires a mirror symmetry in the spanwise direction. The mirror-symmetric solution obtained for $Ra=0$ is identical to the plane Couette flow solutions found in Gibson, Halcrow & Cvitanović (Reference Gibson, Halcrow and Cvitanović2009), Itano & Generalis (Reference Itano and Generalis2009): the mirror-symmetric solution branch for $Re=5000$ is shown by the thin blue dotted curve.

The results shown in figure 1 do not change much for higher Reynolds numbers. For example, the results for $Re=10\,000$ shown by the points almost overlap with the solid and dashed curves. It should be noted that the collapse of the two-dimensional solution branch is not a surprising result, because it is independent of  $Re$ , as remarked earlier. However, in general, the three-dimensional solution does depend on  $Re$ , and the asymptotically convergent result seen in the figure can be regarded as a numerical confirmation of the scaling given in Eaves & Caulfield (Reference Eaves and Caulfield2015). One may notice that the asymptotic convergence of the solutions near the bifurcation point where the mirror-symmetric solution appears is less good. This is because the critical layer becomes nearly flat ( $y=0$ ) as the bifurcation point is approached, and in this special case the asymptotic structure is modified into a localised structure locally described by the boundary region equations (BRE) as concerned in Deguchi & Hall (Reference Deguchi and Hall2015). In fact, when the mirror-symmetric solution is continued to smaller $Ra$ keeping the computational box size, the flow field naturally develops small-scale vortices localised in both $y{-}z$ directions, as the asymptotic theory predicts. The readers who are interested in this regime should refer to Olvera & Kerswell (Reference Olvera and Kerswell2017), where the similar asymptotic development of the solutions is independently found and more carefully examined numerically. As seen in those observations, small-scale vortices are preferred for larger $Ri$ . In the next section, we consider very large wavenumbers comparable to the intrinsic scale of the smallest vortices, in order to preclude the possibility of localisation in  $z$ . In the limit of large wavenumbers it is shown by Deguchi (Reference Deguchi2015) that both of the VWI and BRE states become the UNS states.

Figure 2. The visualisation of the three-dimensional solution indicated by the arrow in figure 1. $(Re,Pr,Ra)=(10\,000,1,-20)$ . (a) Yellow/blue is 50 % maximum/minimum isosurface of streamwise vorticity. (b,c) Streamwise averaged field: (b) is the contour of streamwise velocity, streak, whilst (c) is the contour of temperature.

The flow field of the solution indicated by the arrow in the bifurcation diagram is shown in figure 2. Figure 2(a) shows the isosurface of the streamwise vorticity, which is a good indicator of the wave component. Figures 2(b) and 2(c) are the streak and temperature fields, respectively. These pictures are almost identical since for $Pr=1$ the governing equations for them take the same form at the asymptotic limit. The curved critical layer seen in figure 2(b) suggests that the asymptotic limit of this state is governed by VWI (recall that the critical layer position is defined by $\overline{u}=0$ ). The location where the wave instability is seen in figure 2(a) is in fact the position of the critical layer, as shown for the isothermal cases by Wang, Gibson & Waleffe (Reference Wang, Gibson and Waleffe2007), Hall & Sherwin (Reference Hall and Sherwin2010).

4 Small-scale structures for moderate Richardson numbers

The disrupted flow in Eaves & Caulfield (Reference Eaves and Caulfield2015) is calculated at $Ri=10^{-2}$ , namely $Ra=-10\,000$ . Therefore the scenario concerned in the last section, valid for $Ra\sim O(1)$ , seems to be inapplicable.

Now let us derive the dependence of the typical scale of the vortex $L$ on  $Ra$ . Throughout this section, we select $Pr=1$ , as in Eaves & Caulfield (Reference Eaves and Caulfield2015), for the sake of simplicity. As found by Blackburn et al. (Reference Blackburn, Hall and Sherwin2013), if we choose the wavenumbers to be $O(L^{-1})$ , the VWI states are driven in a layer of vertical thickness $O(L)$ . For given $L$ , we can consider the balance of the buoyancy term in the vortex part as we did in the last section. In the vertical momentum equation in the $x$ -averaged part of (2.3a,b ), the viscous–convective balance requires $O(v)\sim O(1/(Re\,L))$ , whilst the convective–buoyancy balance yields $O(v^{2}/L)\sim O(Ra\,Re^{-2}\unicode[STIX]{x1D703})$ . The final key to the scaling is that outside of the layer where the nonlinear vortices are driven, the flow is predominantly the linear basic flow. The matching of the flows in the inner layer and the outer region requires, as we shall see shortly, the quantities $u/L$ and $\unicode[STIX]{x1D703}/L$ to be $O(1)$ values. Eliminating $v$ and $\unicode[STIX]{x1D703}$ from the viscous–convective balance, the convective–buoyancy balance and $\unicode[STIX]{x1D703}\sim L$ , we arrive at $L\sim O(Ra^{-1/4})$ .

The length scale represents the maximum size of vortices for a given size of  $Ra$ . Thus the size can be related to Ozmidov scale $l_{o}=(\unicode[STIX]{x1D6F7}N^{-3})^{1/2}$ (Ozmidov Reference Ozmidov1965), where the viscous dissipation rate $\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D708}N^{2}$ is estimated by the Brunt–Väisälä frequency  $N$ . In fact, if we note that the Richardson number can be described as $Ri=(N/S)^{2}$ using the dimensional shear $S=u_{0}/h$ , the Ozmidov scale can be expressed as $l_{o}=Ra^{-1/4}h$ . At this scale, the effective Richardson number in the roll equations becomes order unity, and thus the viscous–buoyancy balance is achieved. The roll equation is not directly forced by the basic flow shear, see (3.2a ), and therefore the buoyancy frequency $N$ is the most appropriate choice in estimating  $\unicode[STIX]{x1D6F7}$ .

The localised VWI at each $L$ is self-similar, as shown by Blackburn et al. (Reference Blackburn, Hall and Sherwin2013). However, the smaller the $L$ , the smaller the effective Reynolds number felt by the wave is. Therefore, when the effective Reynolds number becomes of order unity, the problem is no longer amenable to the VWI approach. In that case all terms in the Navier–Stokes equations must be retained, as shown by Deguchi (Reference Deguchi2015). The length scale of this minimal structure is inversely proportional to the square root of the Reynolds number. In fact, it has been shown that this scale corresponds to the Kolmogorov scale $l_{k}=(\unicode[STIX]{x1D6F7}^{-1}\unicode[STIX]{x1D708}^{3})^{1/4}=Re^{-1/2}h$ when the dissipation is estimated by the shear as $\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D708}S^{2}$ . This is the scale at which the effective Reynolds number becomes order unity in the wave equations, and thus the viscous–convective balance is recovered. Here the effect of stratification is absent and thus $S$ is used to estimate $\unicode[STIX]{x1D6F7}$ . If the stratification influences the small scale structure we need to consider the viscous–buoyancy balance as well. In this case we can expect $l_{o}\sim l_{k}$ , from which we arrive at $Ra\sim Re^{2}$ (or equivalently $Ri\sim O(1)$ ).

Now let us rescale the spatial variables as $(x,y,z)=L(X,Y,Z)$ . It is easy to see that if we write

then (2.3) becomes the unit Reynolds number Navier–Stokes equations

where $\unicode[STIX]{x1D735}=(\unicode[STIX]{x2202}_{X},\unicode[STIX]{x2202}_{Y},\unicode[STIX]{x2202}_{Z})$ and

The factor $B/D$ in (4.1), where $B$ is the Nusselt number (i.e. the $x$ – $z$ average of $\unicode[STIX]{x1D703}_{y}|_{y=1}$ ), has been introduced for later convenience. Note that the flow is periodic in $(X,Z)\in [0,2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FC}_{0}]\times [0,2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FD}_{0}]$ with the normalised wavenumbers $\unicode[STIX]{x1D6FC}_{0}=\unicode[STIX]{x1D6FC}L$ and $\unicode[STIX]{x1D6FD}_{0}=\unicode[STIX]{x1D6FD}L$ . In the vertical direction we assume that the all nonlinear interaction is confined to $Y\sim O(1)$ .

Figure 3. Sketch of asymptotic regions for UNS states. The red thick curve is a typical mean flow profile of the corresponding nonlinear solutions. In the inner region with a thickness of $O(Re^{-1/2})$ , the nonlinearity creates a velocity displacement of size $2(D-1)$ for the mean flow. The thin red line is the basic flow. The behaviour of the mean temperature is similar.

The large $|Y|$ behaviour of the inner flow must be chosen to match the outer flow. Since now we assume there is no nonlinearity in the outer regions, the flow there is simply a linear function, see figure 3. Thanks to the factor $B/D$ in (4.1), we can assume the limiting behaviour of the inner flow takes the form

The displacement of velocity and temperature, $d_{0}/\unicode[STIX]{x1D6FD}_{0}$ and $b_{0}/\unicode[STIX]{x1D6FD}_{0}$ respectively, are determined by the solution for given $\unicode[STIX]{x1D6FC}_{0},\unicode[STIX]{x1D6FD}_{0}$ and $R_{0}$ . The factor of $1/\unicode[STIX]{x1D6FD}_{0}$ in the displacements has been introduced so that the definition of $d_{0}$ becomes identical to that used in Deguchi (Reference Deguchi2015).

In the outer region, the streamwise velocity and temperature take the form

where the plus/minus sign must be chosen for the upper/lower outer flow. On the other hand, from the limiting form of the inner flow (4.4), the leading-order outer flow must be

Therefore in order to match (4.5) and (4.6), we require

from which we arrive at a length scale comparable to the Kolmogorov scale discussed earlier. The displacement velocity generated by the nonlinear effect modifies the mean flow as shown in figure 3; the displacement can now be expressed as $2(D-1)=2Re^{-1/2}d_{0}/\unicode[STIX]{x1D6FD}_{0}$ . From (4.7), $R_{0}=Ra\,Re^{-2}(1+\cdots )$ , and thus indeed the small-scale structures persist for large negative Rayleigh numbers of $O(Re^{2})$ , namely $Ri\sim O(1)$ .

This result suggests that the Richardson number at $Ri_{c}(Re)$ may tend to a constant for large Reynolds number. One may notice that this conclusion differs from the previously known results, where it was argued that strongly stratified turbulent flows are governed by a buoyancy Reynolds number ${\sim}O(Re/Ri)$ . If such asymptotic states exist, $Ri_{c}(Re)$ should ultimately be proportional to $Re$ ; see Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), for example. However, their numerical simulations are subjected to certain random forcing to inject energy to the system, and thus the flow configuration is not exactly the same as the one considered in this paper.

Figure 4. The bifurcation diagram of large wavenumber plane Couette flow solutions for $Pr=1$ , $\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D6FC}=2$ . The normalised displacement velocity $d_{0}$ is shown against $\unicode[STIX]{x1D6FD}_{0}^{2}/4$ (that quantity should have appeared in the horizontal axis of figure 2(c) in Deguchi (Reference Deguchi2015)). Here the scaled spanwise wavenumber $\unicode[STIX]{x1D6FD}_{0}\approx Re^{-1/2}\unicode[STIX]{x1D6FD}$ . Curves and points are computed at $Re=5000$ and $Re=10\,000$ , respectively. Dashed curve: isothermal case $Ra=0$ . Solid curve: stably stratified case $Ra=-0.001Re^{2}$ .

In figure 4 we confirm the above scaling in the large wavenumber computation of the steady solution. For the fixed wavenumber ratio $\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D6FC}=2$ , the value of the normalised displacement $d_{0}$ is monitored while changing the scaled wavenumber $\unicode[STIX]{x1D6FD}_{0}$ . Here, the curves and points are computed for $Re=5000$ and $Re=10\,000$ , respectively. The dashed curve and the open circles are computed for the isothermal case $Ra=0$ ; the results are essentially those given in figure 2(c) of Deguchi (Reference Deguchi2015). Similar excellent asymptotic convergence can be found for the stably stratified case $Ra=-0.001Re^{2}$ , shown by the solid curve and the filled circles. The flow structure of the solution at $(Re,Ra)=(10\,000,-100\,000)$ , corresponding to the solution indicated by the arrow in figure 4, is shown in figure 5. The streamwise vorticity, streamwise velocity perturbation, and temperature perturbation are localised near the centre of the channel, consistent with the theory described here. In our computation, the perturbation vortices should appear at the centre of the channel due to the symmetry assumed. If there is no symmetry restriction, we may be able to find localised structures for other vertical positions, but such structures are merely a Galilean transformed version of the localised structure seen in figure 5. In fact, similar to the isothermal cases, our localised solution can be used as a good asymptotic approximation of locally appeared coherent structures in any stratified shear flow, as long as we can assume the background velocity and temperature profiles can be locally approximated by a linear profile.

Figure 5. The same plot as figure 2 but for the solution indicated by the arrow in figure 4. $(Re,Pr,Ra,\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD})=(10\,000,1,-100\,000,6.3,12.6)$ .