6.1 Shallow horizontal layer $(A\gg 1)$

In this subsection, the effects of the governing parameters on the convective flows, as predicted by the present analytical solution, are compared with the numerical results obtained by solving the full governing equations. As discussed above, the governing parameters of the problem are the Rayleigh number $R_{T}$ , the buoyancy ratio $\unicode[STIX]{x1D711}$ , the Lewis number $Le$ , the normalized porosity $\unicode[STIX]{x1D700}$ and the modified Prandtl number $P_{r}^{\ast }$ .

Figure 4 exemplifies the effect of the modified Prandtl number, $P_{r}^{\ast }$ , on the onset of motion and the resulting convective flow in the present system. The form drag effect is proportional to $1/P_{r}^{\ast }$ . The bifurcation diagrams are presented in terms of the flow intensity, $\unicode[STIX]{x1D6F9}_{0}$ , the value of the stream function at the centre of the cavity, as a function of $R_{T}$ for $\unicode[STIX]{x1D711}=-0.1$ , $Le=5$ and $\unicode[STIX]{x1D700}=1$ . The curves depicted in the graphs are the predictions of the present analytical/numerical nonlinear theories. The solid lines correspond to stable branches and the dot-dot-dashed lines to unstable ones. The numerical solutions of the full governing equations, obtained for $A=10$ , are shown by solid circles. A good agreement is observed between these two nonlinear theory results. For the values of the governing parameters considered here, the critical Rayleigh numbers for the onset of motion, as predicted by (4.22), are given by $R_{TC}^{sup}=24$ , $R_{TC}^{over}=16$ and $R_{TC}^{osc,0}=20.53$ . These critical Rayleigh numbers are not affected by the form drag since they correspond to onsets of motion occurring from the rest state. It is well known that the flow circulation in the layer can be indifferently clockwise or counterclockwise. Since the results are symmetrical, only the counterclockwise solutions $(\unicode[STIX]{x1D6F9}_{0}>0)$ are depicted in all of the bifurcation curves presented in this paper.

The results presented in figure 4(a) are obtained for $P_{r}^{\ast }\rightarrow \infty$ , i.e. in the absence of form drag effects. The resulting bifurcation curve is consistent with the behaviour reported in the past by Mamou & Vasseur (Reference Mamou and Vasseur1999) and Bahloul et al. (Reference Bahloul, Vasseur and Robillard2007), while considering double-diffusive and Soret-induced convection in a porous layer modelled according to Darcy’s law. The onset of motion is characterized by a transition from the rest state through finite amplitude convection $(\unicode[STIX]{x1D6F9}_{0}\neq 0)$ , i.e. nonlinear instabilities. This transition occurs at a subcritical Rayleigh number of $R_{TC}^{sub}=15.63$ , which is well below the supercritical Rayleigh number of $R_{TC}^{sup}=24$ . At this threshold, subcritical bifurcation occurs only when the buoyancy forces induced by the thermal and solutal contributions are opposing each other and when the stabilizing agent is the slower diffusing component $(Le>1)$ . In the graph, the solution depicted by the solid line (thermal branch), referred to as the stable solution, is confirmed by the numerical results depicted by black circles. This is not the case for the dot-dot-dashed line (solutal branch), connecting the subcritical and supercritical Rayleigh numbers, which corresponds to an unstable solution. The thermal and the solutal lines are connected to each other at a saddle–node bifurcation. In general, the numerical results presented in figure 4(a) indicate that, below the subcritical Rayleigh number (region I), the rest state prevails. In region (II), delineated by $R_{TC}^{sub}\leqslant R_{T}\leqslant R_{TC}^{over}$ , parallel flow theory predicts the existence of two solutions (one stable and the other one not). In region (III), delimited by $R_{TC}^{over}\leqslant R_{T}\leqslant R_{TC}^{osc,0}$ , starting with the conduction state as initial conditions, the numerical results indicate that for $R_{T}=16.1$ the flow is oscillatory. The result is presented in the graph, by circle symbols, in terms of the root mean square of the stream function, $\overline{\unicode[STIX]{x1D6F9}}_{0}$ . Upon increasing the Rayleigh number step by step and using the previous results as initial conditions, the oscillating flow branch $\overline{\unicode[STIX]{x1D6F9}}_{0}$ can be sustained up to $R_{T}=16.9$ . Above this Rayleigh number, the solution bifurcates directly towards the steady convective motion predicted by parallel flow theory. Numerically, it was found that upon increasing (decreasing) the Rayleigh number and starting the computations with a conductive state (finite amplitude convection) the resulting solution follows the hysteresis loop depicted by the arrow paths depicted on the graphs.

The results presented in figure 4(b) are obtained for a relatively small value of the modified Prandtl number, namely $P_{r}^{\ast }=1$ , for which the form drag effect is significant. As already mentioned, the critical Rayleigh numbers for the onset of motion from the rest state ( $R_{TC}^{over}$ , $R_{TC}^{osc,0}$ , $R_{TC}^{sup}$ ) are not affected by the form drag effect. Thus, the minimum critical Rayleigh number for the onset of linear motion is the overstable Rayleigh number $R_{TC}^{over}=16$ , such that below this value, zone (I), the system remains at rest, as indicated by the numerical results presented in the graph. Above this value, similarly to the results obtained in figure 4(a), for $R_{T}>R_{TC}^{over}$ the numerical results reveal the existence of an oscillating flow branch, $\overline{\unicode[STIX]{x1D6F9}}_{0}$ , over the range $16.2\leqslant R_{T}\leqslant 28.5$ . The resulting type of oscillating flow, in zone (II), is illustrated in figure 4(b) for the particular case $R_{T}=18$ in terms of the time history of $\unicode[STIX]{x1D6F9}_{0}$ for a few cycles. It is found that the oscillating branch also exists in zone (III), delineated by $R_{TC}^{osc,0}\leqslant R_{T}\leqslant R_{TC}^{sup}$ , as exemplified by the results obtained for $R_{T}=22$ . The oscillatory solution branch is also observed to extend beyond $R_{T}=R_{TC}^{sup}$ , where it bifurcates toward the steady state solution branch. The solid line in zone (IV) is the prediction of nonlinear parallel flow theory. It is observed that the onset of steady motion is now not subcritical, compared with figure 4(a) for Darcy’s law, but rather supercritical. The resulting pitchfork type bifurcation, characterized by a transition from the rest state $(\unicode[STIX]{x1D6F9}_{0}=0)$ to a convective regime, takes place at a supercritical Rayleigh number of $R_{TC}^{sup}=24$ . Upon increasing the Rayleigh number above this critical value, a good agreement is observed between the analytical and numerical results.

Figure 5.

Two possible time histories of the stream function intensity $\unicode[STIX]{x1D6F9}_{0}$ and the flow structure for $\unicode[STIX]{x1D711}=-0.1$ , $Le=5$ , $P_{r}^{\ast }=2$ and $\unicode[STIX]{x1D700}=1$ : (a)  $R_{T}=18$ (zone II), (b)  $R_{T}=22$ (zone III).

The case $P_{r}^{\ast }=2$ , which is an intermediate value of the modified Prandtl number between those considered above, is depicted in figure 4(c). The bifurcation diagram is now characterized by six different zones. In particular, zones (I) and (II), delineated by $R_{TC}^{over}$ and $R_{TC}^{osc,0}$ , are the same as those discussed above, for $P_{r}^{\ast }=1$ . Thus, the rest state prevails in zone (I), while for $R_{T}>R_{TC}^{over}$ the flow is found to be oscillatory. This flow regime gives rise to a branch corresponding to the numerical results obtained for $\overline{\unicode[STIX]{x1D6F9}}_{0}$ versus $R_{T}$ . It is found that the oscillating branch can be maintained up to $R_{T}=23.9$ , above which the numerical solution bifurcates towards a steady finite amplitude convective regime. Figure 5(a) illustrates the results obtained for $R_{T}=18$ , zone (II), in terms of the variation of the stream function $\unicode[STIX]{x1D6F9}_{0}$ in time and selected streamline patterns during a period. The flow structure consists of a unicellular cell filling up the entire cavity. An overturning of the flow from counterclockwise to clockwise circulation and vice versa is observed. The bifurcation curve, predicted by parallel flow theory, indicates the onset of motion occurring now at a subcritical Rayleigh number of $R_{TC}^{sub}=23.51$ . Thus, zone (III) in the graph is now confined in the range $R_{TC}^{osc,0}\leqslant R_{T}\leqslant R_{TC}^{sub}$ . Figure 5(b) illustrates the results obtained in this region, for $R_{T}=22$ . The strength of the resulting oscillating circulation is found to be naturally higher than that reported in zone (II) and the corresponding flow structure is more complex. At time $t=620$ , the flow pattern consists of a counterclockwise unicellular circulation. As time increases, this unicellular structure breaks up to give rise to a central clockwise circulation and two small counterclockwise rotating eddies, localized near the end regions of the cavity. These two eddies vanish gradually, giving rise to a weak unicellular clockwise circulation in the layer, as can be seen from the results obtained at $t=840$ . As time progresses further, the flow pattern consists of a counterclockwise circulating cell in the central part of the system and two clockwise eddies near the end regions. At time $t=1020$ , the unicellular counterclockwise circulation observed at $t=620$ is recovered. The whole process is repeated during each time period. This oscillatory behaviour and the flow patterns could be different when the layer aspect ratio is varied. The results, predicted by the present parallel flow model, are now discussed. The resulting bifurcation curve is seen to be quite particular. A bird’s-eye view of the graph indicates the existence of two stable different branches, depicted by solid lines. The upper one, starting at $R_{TC}^{sub}$ , is the standard branch observed in subcritical bifurcations, as exemplified, for instance, in figure 4(a). On the other hand, the lower stable branch, which starts at $R_{TC}^{sup}=24$ , corresponds to a supercritical bifurcation similar to that discussed in figure 4(b). However, it is observed that this latter exists only in the range $R_{TC}^{sup}\leqslant R_{T}\leqslant R_{TC}^{tur}$ , where $R_{TC}^{tur}=25.08$ corresponds to a turning saddle–node point. The two stable branches are connected with each other by an unstable branch, shown as a dot-dot-dashed line. This is a typical bistability bifurcation. Numerical confirmation of the analytical solution is also presented in the graph and here again the agreement is excellent. Starting the numerical procedure with the rest state as initial conditions, it was found that, upon increasing the Rayleigh number above $R_{TC}^{sub}=23.51$ , the resulting solution follows the hysteresis loop depicted by arrows in the zoom included in figure 4(c). Thus, for $R_{T}\geqslant R_{TC}^{tur}$ , zone (VI), only the upper branch exists and it is seen that, as expected, upon increasing $R_{T}$ the strength of the convection is promoted monotonically. Such a steady solution can be maintained up to a critical Rayleigh number of $R_{TC}^{Hopf}$ , above which a Hopf bifurcation occurs and the flow becomes unsteady. Hopf bifurcation will be covered in the following section. On the other hand, upon decreasing the Rayleigh number below $R_{TC}^{tur}$ , one must now consider two distinct zones. Zone (V) is characterized by the existence for a given Rayleigh number of three possible different solutions (two stable and one unstable). Region (V), which is demarcated by $R_{TC}^{sup}\leqslant R_{T}\leqslant R_{TC}^{tur}$ , is the cradle where a bistability phenomenon occurs leading to the existence of two distinct finite amplitude solutions between two forward and backward subcritical bifurcations, even though it happens beyond the thresholds of supercritical convection. Within this region, two lower and upper stable convective branches were obtained and connected by a third unstable branch.

Figure 6(a,b) illustrates the streamlines, isotherms and isoconcentrates obtained numerically in the bistability convection zone, for $R_{T}=24.5$ . Appropriate initial conditions had to be specified to simulate the two different stable solutions. Thus, to help in bifurcating towards the desired solution, the analytical solution was used as initial condition. For the upper branch, the lower branch can be obtained naturally when departing from the rest state. The upper-branch flow patterns are shown in figure 6(a) and the lower-branch patterns are displayed in figure 6(b). Naturally, it has not been possible to confirm numerically the third unstable solution predicted by the analytical model.

Figure 6.

Contours of stream (top), temperature (middle) and concentration (bottom) for $R_{T}=24.5$ , $P_{r}^{\ast }=2$ , $\unicode[STIX]{x1D711}=-0.1$ , $Le=10$ and $\unicode[STIX]{x1D700}=1$ ; (a)  $\unicode[STIX]{x1D6F9}_{0}=0.394$ , $Nu=1.063$ , $Sh=2.182$ ; (b)  $\unicode[STIX]{x1D6F9}_{0}=0.016$ , $Nu=1.000$ , $Sh=1.003$ .

Figure 7.

Bifurcation diagram as a function of $R_{T}$ and $P_{r}^{\ast }$ for $\unicode[STIX]{x1D711}=-0.1$ , $Le=5$ and $\unicode[STIX]{x1D700}=1$ ; (a) flow intensity $\unicode[STIX]{x1D6F9}_{0}$ , (b) Nusselt number $Nu$ and (c) Sherwood number $Sh$ ; ——, analytical stable;

, analytical unstable;

, numerical.

A more complete view of the effects of $P_{r}^{\ast }$ and $R_{T}$ on $\unicode[STIX]{x1D6F9}_{0}$ , $Nu$ and $Sh$ is presented in figure 7(a–c) for $\unicode[STIX]{x1D711}=-0.1$ , $Le=5$ and $\unicode[STIX]{x1D700}=1$ . Results are obtained for values of the modified Prandtl number ranging between $P_{r}^{\ast }=\infty$ (Darcy) and 1. In general, it is seen from figure 7(a) that, for a given Rayleigh number, a decrease in the form drag parameter reduces the strength of the convective motion in the layer considerably. This is due to the fact that decrease of $P_{r}^{\ast }$ results in an increase of the form drag. Consequently, as can be observed from figure 7(b,c), both the heat and the mass transfer rates are considerably inhibited with a decrease of the modified Prandtl number. The plots also clearly show the effect of $P_{r}^{\ast }$ on the existence of the different types of bifurcation curves discussed in figure 4. For the set of governing parameters considered here, it is predicted by the analytical model that the bifurcation is subcritical provided that $P_{r}^{\ast }>1.85$ , as depicted in figure 8(a) for $P_{r}^{\ast }=2$ . The resulting bifurcation curve is characterized by an upper subcritical stable branch existing for $R_{T}\geqslant R_{TC}^{sub}=23.5$ . Furthermore, due to the form drag effect, the graph also indicates the existence of a supercritical stable branch extending in the range $R_{TC}^{sup}=24\leqslant R_{T}\leqslant R_{TC}^{tur}=25.1$ . Both $R_{TC}^{sub}$ and $R_{TC}^{tur}$ are connected to each other by an unstable branch depicted by a dot-dot-dashed line. Upon decreasing the value of $P_{r}^{\ast }$ , the onset of motion remains subcritical down to $P_{r}^{\ast }=1.85$ , at which point the condition $R_{TC}^{sub}=R_{TC}^{sup}=24$ is reached, as illustrated in figure 8(b). For $P_{r}^{\ast }\leqslant 1.85$ , it is seen from figure 8(b,c) that the onset of motion is now supercritical ( $R_{TC}^{sub}\geqslant R_{TC}^{sup}$ ). Thus, as exemplified in figure 8(c) for $P_{r}^{\ast }=1.5$ , it is noticed that the subcritical Rayleigh number $R_{TC}^{sub}$ becomes greater than $R_{TC}^{sup}$ , such that it degenerates into a second turning point $R_{TC}^{sub}$ . The bifurcation curve is now characterized by two turning points, $R_{TC}^{sub}=25.4$ and $R_{TC}^{tur}=25.9$ . It is noted from figure 8(a–c) that, in the range $R_{TC}^{sup}(R_{TC}^{sub})\leqslant R_{T}\leqslant R_{TC}^{tur}$ , three different solutions (two stable and one unstable) are possible for a given $R_{T}$ , confirming the existence of the bistability phenomenon discussed above. Finally, the results presented in figure 8(d) for $P_{r}^{\ast }=1.3$ indicate that the condition $R_{TC}^{tur}=R_{TC}^{sub}=26.5$ is reached. As a result, for $P_{r}^{\ast }\leqslant 1.3$ , the analytical model predicts that, for a given Rayleigh number $R_{T}>R_{TC}^{sup}$ , the resulting supercritical bifurcation curves are now single valued. A stability diagram is illustrated in figure 8(e), where the bistability region is delineated by the hatched area.

Figure 8.

Bifurcation diagram in terms of $\unicode[STIX]{x1D6F9}_{0}$ versus $R_{T}$ for $\unicode[STIX]{x1D711}=-0.1$ and $Le=5$ ; (a)  $P_{r}^{\ast }=2$ , (b)  $P_{r}^{\ast }=1.85$ , (c)  $P_{r}^{\ast }=1.5$ and (d)  $P_{r}^{\ast }=1.3$ ; ——, analytical stable;

, analytical unstable. (e) Stability diagram showing the bistability region (hatched) and the position of the bifurcations related to graphs (a–d).

Figure 9.

Analytical results obtained for $P_{r}^{\ast }=2$ , $Le=5$ and $\unicode[STIX]{x1D700}=1$ : (a) bifurcation diagram in terms of $\unicode[STIX]{x1D6F9}_{0}$ versus $R_{T}$ , (b) critical Rayleigh number $R_{TC}$ versus  $\unicode[STIX]{x1D711}$ .

Figure 9(a) shows typical bifurcation diagrams for $\unicode[STIX]{x1D6F9}_{0}$ as a function of $R_{T}$ and $\unicode[STIX]{x1D711}$ for $P_{r}^{\ast }=2$ , $Le=5$ and $\unicode[STIX]{x1D700}=1$ . The curve obtained for $\unicode[STIX]{x1D711}=0$ corresponds to convection induced solely by the thermal effect. It gives rise to a pitchfork bifurcation occurring at a supercritical Rayleigh number of $R_{TC}^{sup}=12$ . The way in which this supercritical bifurcation is affected upon decreasing the value of the buoyancy ratio is illustrated in the graph. In fact, the evolution of the bifurcation curves with the buoyancy ratio is quite similar to that discussed in figure 8(a–d) while studying the effect of $P_{r}^{\ast }$ . Figure 9(b) shows the variation of the critical Rayleigh numbers $R_{TC}^{sup}$ , $R_{TC}^{sub}$ and $R_{TC}^{tur}$ with $\unicode[STIX]{x1D711}$ . The graph indicates that in the range $-0.15\leqslant \unicode[STIX]{x1D711}\leqslant -0.096$ , the subcritical Rayleigh number is the lower critical Rayleigh number $(R_{TC}^{sub}<R_{TC}^{sup}<R_{TC}^{tur})$ , such that the onset of motion is subcritical, as exemplified in figure 8(a). At $\unicode[STIX]{x1D711}=-0.096$ , it is found that $R_{TC}^{sub}=R_{TC}^{sup}=23.19$ , this situation being similar to that depicted in figure 8(b). Upon keeping on decreasing the value of the buoyancy ratio, it is observed that, in the range $-0.096\leqslant \unicode[STIX]{x1D711}\leqslant -0.081$ , the bifurcation curve consists of a lower supercritical branch extending in the range $R_{TC}^{sup}\leqslant R_{T}\leqslant R_{TC}^{tur}$ and an upper branch that exists for the condition $R_{T}\geqslant R_{TC}^{sub}$ , as exemplified in figure 8(c). The subcritical branch has thus degenerated into a branch corresponding to the critical Rayleigh number $R_{TC}^{sub}$ . The results obtained at $\unicode[STIX]{x1D711}=-0.081$ , for which $R_{TC}^{tur}=R_{TC}^{sub}=21.58$ , are similar to those depicted in figure 8(d), i.e. the two turning point branches merge. For $-0.081\leqslant \unicode[STIX]{x1D711}\leqslant 0$ , the only existing critical Rayleigh number is the supercritical one, $R_{TC}^{sup}$ .

Figure 10.

(a) Subcritical Rayleigh number $R_{TC}^{sub}$ as a function of $P_{r}^{\ast }$ and $Le$ for $\unicode[STIX]{x1D711}=-0.1$ . (b) Bifurcation diagram in terms of $\unicode[STIX]{x1D6F9}_{0}$ as a function of $R_{T}$ and $P_{r}^{\ast }$ for $\unicode[STIX]{x1D711}=-0.1$ , $Le=15$ and $\unicode[STIX]{x1D700}=1$ ; ——, analytical stable;

, analytical unstable;

, numerical.

The subcritical Rayleigh number $R_{TC}^{sub}$ is plotted in figure 10(a) as a function of $P_{r}^{\ast }$ and various values of the Lewis number, for the case $\unicode[STIX]{x1D711}=-0.1$ . For larger values of the modified Prandtl number $(P_{r}^{\ast }\rightarrow \infty )$ , the results predicted by Darcy’s model, equation (4.23), are recovered, as depicted in the small inset included in the graph. Thus, as $Le\rightarrow \infty$ , the solute concentration is almost uniform in the system, except in very thin layers near the horizontal boundaries. Consequently, the onset of convection approaches asymptotically the criterion for the onset of convection in a porous layer subjected solely to constant heat fluxes (i.e.  $R_{TC}^{sub}\rightarrow R_{TC}^{sup}=12$ ). Upon decreasing the Lewis number progressively, $R_{TC}^{sub}$ first increases sharply, up to a maximum, $R_{TC}^{sub}=16$ at $Le=3.22$ , and then decreases monotonically. It is observed that below $Le=2.16$ , subcritical convection is not possible. This follows from the fact that when the Lewis number is small enough, the mass transfer process across the porous layer is ruled by pure diffusion. The onset of motion is then supercritical, occurring at a Rayleigh number of $R_{TC}^{sup}=15.3$ , as predicted by (4.22). Upon keeping on decreasing the Lewis number further, it is seen that the supercritical Rayleigh number tends asymptotically towards $R_{TC}^{sup}=12$ as $Le\rightarrow 0$ . The effect of the modified Prandtl number on the subcritical Rayleigh number is clearly illustrated in the graph. As already mentioned, a decrease of $P_{r}^{\ast }$ annihilates the strength of the convective motion in the layer. Thus, for a given Lewis number, it is found that the value of $R_{TC}^{sub}$ , at which convection occurs, increases considerably upon decreasing $P_{r}^{\ast }$ . For instance, for $Le=10^{3}$ , it is found that the subcritical Rayleigh number increases from $R_{TC}^{sub}=12.25$ for $P_{r}^{\ast }=10^{3}$ to $R_{TC}^{sub}=28$ for $P_{r}^{\ast }=10^{-1}$ . Naturally, for the lower values of the Lewis number presented in the graph, subcritical convection is not always possible as the modified Prandtl number is decreased. Hence, for $Le=3(5)$ , the onset of motion is supercritical for $P_{r}^{\ast }\leqslant 14(1.85)$ . Another view of the effect of the modified Prandtl number is illustrated in figure 10(b). The results obtained in terms of $\unicode[STIX]{x1D6F9}_{0}$ as a function of $R_{T}$ and $P_{r}^{\ast }$ , for $\unicode[STIX]{x1D711}=-0.1$ and $Le=15$ , are found to be different from those discussed in figure 7(a) for $Le=5$ . Thus, independently of the value of $P_{r}^{\ast }$ , the analytical solution indicates that the subcritical bifurcation curves consist of upper stable branches connected with lower unstable branches. The latter are seen to tend asymptotically towards the supercritical Rayleigh number $R_{TC}^{sup}\rightarrow \infty$ . In fact, for the particular Lewis number considered here, onset of supercritical convection is possible only for a negative value of the Rayleigh number, namely $R_{TC}^{sup}=-24$ , as predicted by (4.22). A negative value of $R_{T}$ corresponds to the case of a layer heated from the top. Decrease of the value of the Rayleigh number below this critical value yields the pitchfork bifurcation curves depicted in the graph. Although the strength of the resulting convective motions $\unicode[STIX]{x1D6F9}_{0}$ is very small, the enlargement depicted in figure 10(b) shows that the agreement between the analytical model and the numerical code is very good. For a given value of the form drag coefficient, it is seen that the rest state prevails in the range $-24\leqslant R_{T}\leqslant R_{TC}^{sub}(P_{r}^{\ast })$ .

Figure 11(a) shows the influence of the form drag coefficient on the thresholds of bifurcations $(R_{TC}^{Hopf},R_{TC}^{sub},R_{TC}^{sup},R_{TC}^{over},R_{TC}^{osc,0})$ for $R_{S}=-240$ , $Le=10$ and $\unicode[STIX]{x1D700}=0.2$ . According to (4.22), the onsets of motion from the rest state are given by $R_{TC}^{over}=138$ , $R_{TC}^{osc,0}=180$ and $R_{TC}^{sup}=252$ , these values being independent of $P_{r}^{\ast }$ . As expected, the graph indicates that the onset of Hopf bifurcation, $R_{TC}^{Hopf}$ , decreases sharply upon increasing the value of the form drag parameter. It follows that the steady parallel flow is destabilized earlier with an increase of $P_{r}^{\ast }$ . A similar trend is observed for the evolution of the subcritical Rayleigh number, $R_{TC}^{sub}$ , with the form drag coefficient. In both cases, for sufficiently large values of $P_{r}^{\ast }$ , the Darcy limits represented by dashed lines are reached asymptotically. It is also observed that, in the range $0.12\leqslant P_{r}^{\ast }\leqslant 10^{3}$ , the minimum critical Rayleigh number for the onset of motion is the subcritical one. However, for $P_{r}^{\ast }\leqslant 0.12$ , linear stability theory predicts that convection is possible when $R_{T}\geqslant R_{TC}^{over}$ , the overstable critical Rayleigh number. Figure 11(b) illustrates the effect of $P_{r}^{\ast }$ on the wavelength $A_{C}$ of the incipient flow structure at the onset of Hopf bifurcation. It is found that the wavelength (frequency) decreases (increases) monotonically from $A_{C}=1.35$ ( $f_{r}=21$ ), in the pure Darcy regime, towards $A_{C}=1.22$ ( $f_{r}=22.67$ ), in the pure Dupuit regime. The form drag effects on these two parameters are thus relatively weak.

Figure 11.

Form drag effect for $R_{S}=-240$ , $Le=10$ , $Pr^{\ast }=10$ , and $\unicode[STIX]{x1D700}=0.2$ on (a) the thresholds of convection $R_{TC}$ and (b) the critical wavelength $A_{C}$ and oscillatory frequency  $f_{r}$ .

Figure 12.

Bifurcation diagram in terms of $\unicode[STIX]{x1D6F9}_{0}$ versus $R_{T}$ for $R_{S}=-240$ , $Le=10$ , $Pr^{\ast }=10$ , and $\unicode[STIX]{x1D700}=0.2$ ; ——, analytical stable;

, analytical unstable; ●, numerical stable;

, numerical unsteady.

Figure 12 exemplifies the effects of the thermal Rayleigh number, $R_{T}$ , on the flow intensity, $\unicode[STIX]{x1D6F9}_{0}$ , for the case where $R_{S}=-240$ , $P_{r}^{\ast }=10$ , $Le=10$ and $\unicode[STIX]{x1D700}=0.2$ . The numerical results obtained for a cavity with an aspect ratio of $A=20$ are observed to be in excellent agreement with the prediction of the parallel flow model. According to (4.22), the critical Rayleigh numbers for the onset of motion from the rest state are given by $R_{TC}^{over}=138$ , $R_{TC}^{osc,0}=180$ and $R_{TC}^{sup}=252$ . The onset of oscillating motion, predicted by the numerical procedure described in § 3, occurs at $R_{TC}^{Hopf}=1967$ , while the onset of subcritical convection predicted by the nonlinear parallel flow approximation is $R_{TC}^{sub}=29$ . In the graph, zones (I)–(V) are equivalent to those defined in figure 4(a), while in zone (VI), for which $R_{T}\geqslant R_{TC}^{Hopf}$ , the convective flow is oscillating. This point is confirmed by the numerical results obtained for $R_{T}=2200$ , 3000 and 4000. The empty circles presented in the graph correspond to the resulting time-averaged values of $\unicode[STIX]{x1D6F9}_{0}$ obtained by the numerical code. The time evolution of the flow intensity, $\unicode[STIX]{x1D6F9}_{0}$ , obtained for $R_{T}=2200$ , just above the critical Rayleigh number for Hopf bifurcation predicted by linear stability theory $(R_{TC}^{Hopf}=1967)$ , is depicted in figure 13(a). The resulting oscillatory flow is found to be simply periodic. Snapshots of the perturbations of the stream function (obtained from subtracting the time-averaged solution from the time-accurate one) at different time steps shown in figure 13(a) are presented in figure 13(b). The numerical results indicate that the critical wavelength and the oscillation frequency are given by $A_{C}=1.25$ and $f_{r}=23.15$ respectively, which are close to the values $A_{C}=1.24$ and $f_{r}=22.30$ predicted by linear stability theory. At the onset of Hopf bifurcation, as predicted by the linear stability analysis, there exist two symmetrical solutions at the same circular frequency $p_{i}$ (imaginary part of $p$ ) but with opposite signs. The negative and positive signs indicate the existence of travelling waves in the horizontal direction in both ways, as demonstrated by the full numerical solution presented in figure 13(b) and by the supplementary movies (available at https://doi.org/10.1017/jfm.2016.787). The incipient perturbation flow patterns, obtained by the linear stability analysis, are illustrated in figure 11(b), and once superposed they appear to be similar to those obtained numerically in figure 13(b).

Figure 13.

(a) Time history of the stream function $\unicode[STIX]{x1D6F9}_{0}$ as predicted by the numerical solution and (b) snapshots of the perturbations of the stream function $\unicode[STIX]{x1D6F9}_{p}$ corresponding to selected points, for $R_{T}=2200$ , $R_{S}=-240$ , $Le=10$ , $Pr^{\ast }=10$ , $A=20$ and $\unicode[STIX]{x1D700}=0.2$ . Animations of the stream function, temperature and solute concentration perturbations can be seen in the supplementary movies.

6.2 Square cavity $(A=1)$

In this section, we consider the case of a square cavity subject to vertical temperature and concentration gradients. The two vertical walls of the cavity are impermeable and adiabatic, while Dirichlet or Neumann boundary conditions on temperature and solute are imposed on the horizontal walls. For this situation, an analytical model is not available, such that all of the results presented here are the predictions of the numerical code discussed in § 3.

Figure 14.

Bifurcation diagram in terms of $\unicode[STIX]{x1D6F9}_{0}$ as a function of $R_{T}$ and $P_{r}^{\ast }$ for $A=1$ , $\unicode[STIX]{x1D711}=-0.1$ , $Le=5$ and $\unicode[STIX]{x1D700}=1$ ; (a)  $P_{r}^{\ast }=\infty$ –1; (b)  $P_{r}^{\ast }=2$ ;

, unstable; 

, numerical. (c) Contours of stream (left), temperature (middle) and concentration (right) for $R_{T}=48$ , $P_{r}^{\ast }=2$ : $\unicode[STIX]{x1D6F9}_{0}=0.433$ , $Nu_{m}=1.040$ , $Sh_{m}=1.779$ (upper solution). (d) Contours of stream (left), temperature (middle) and concentration (right) for $R_{T}=48$ , $P_{r}^{\ast }=2$ : $\unicode[STIX]{x1D6F9}_{0}=0.047$ , $Nu_{m}=1.000$ , $Sh_{m}=1.013$ (lower solution). (Neumann boundary conditions.)

Figure 14(a–d) shows typical results obtained for the case of a square cavity heated and salted from the bottom by constant fluxes. Figure 14(a) illustrates the bifurcation curve, depicting the effects of the Rayleigh number, $R_{T}$ , and the modified Prandtl number, $P_{r}^{\ast }$ , on the flow intensity, $\unicode[STIX]{x1D6F9}_{0}$ , for the case $\unicode[STIX]{x1D711}=-0.1$ , $Le=5$ and $\unicode[STIX]{x1D700}=1$ . The solid lines in the graph are obtained by a curve fitting of the numerical data, while the dot-dot-dashed lines are manually guessed predictions (artistic) based on the results obtained for the shallow layer. A bird’s-eye view of the graph indicates that the results are quantitatively similar to those obtained in the previous section for a high-aspect-ratio layer. Thus, upon decreasing the value of the Prandtl number from $P_{r}^{\ast }=\infty$ to 1, it is seen that the bifurcation curves evolve from a subcritical behaviour to a supercritical one. For intermediates values of the modified Prandtl number, the bifurcation curves indicate the occurrence of stable supercritical and subcritical branches connected by an unstable one. This point is exemplified by the curve for $P_{r}^{\ast }=2$ depicted in figure 14(b), showing again the region where bistability convection occurs. For this situation, two finite and distinct convective solutions are possible for a given range of the Rayleigh number depending upon the governing parameters of the problem. This point is illustrated in figure 14(c,d), showing the streamlines, isotherms and isoconcentrates obtained for $R_{T}=48$ . The solution presented in figure 14(d), which corresponds to the lower branch, is obtained starting the numerical calculations using the rest state as initial conditions. However, the solution corresponding to the upper branch, figure 14(c), can be obtained only upon following the hysteresis path depicted on the curve for $P_{r}^{\ast }=2$ . Figure 15 shows typical results obtained for a cavity subject to Dirichlet type thermal and solutal boundary conditions. The governing parameters considered here are the same as those for figure 14. It is clear that the existence of dual solutions is also obtained, as depicted by the bifurcation curve obtained, for instance, for $P_{r}^{\ast }=5$ .

Figure 15.

Bifurcation diagram in terms of $\unicode[STIX]{x1D6F9}_{0}$ as a function of $R_{T}$ and $P_{r}^{\ast }$ for $A=1$ , $\unicode[STIX]{x1D711}=-0.1$ , $Le=5$ and $\unicode[STIX]{x1D700}=1$ ; (a)  $P_{r}^{\ast }=\infty$ –1; (b)  $P_{r}^{\ast }=5$ ;

, unstable; 

, numerical. (c) Contours of stream (left), temperature (middle) and concentration (right) for $R_{T}=80$ , $P_{r}^{\ast }=5$ : $\unicode[STIX]{x1D6F9}_{0}=1.005$ , $Nu_{m}=1.108$ , $Sh_{m}=2.274$ (upper solution). (d) Contours of stream (left), temperature (middle) and concentration (right) for $R_{T}=80$ , $P_{r}^{\ast }=5$ : $\unicode[STIX]{x1D6F9}_{0}=0.028$ , $Nu_{m}=1.000$ , $Sh_{m}=1.002$ (lower solution). (Dirichlet boundary conditions.)

The multiplicity of solutions and bistability convection reported in this study is due to the effects of the porous medium form drag $(P_{r}^{\ast })$ on the subcritical bifurcation flow regimes that occur when the buoyancy forces induced by the thermal and solutal contributions are opposing each other $(\unicode[STIX]{x1D711}<0)$ and when the stabilizing agent is the slower diffusing component $(Le>1)$ . Due to the symmetry of the boundary conditions considered here, the circulation of the convective flow in the cavity rotates indifferently clockwise or counterclockwise. All of the bifurcation curves presented in this paper are for the counterclockwise solutions $(\unicode[STIX]{x1D6F9}_{0}>0)$ , as exemplified, for instance, by figure 4(c). It follows that, for appropriate values of the governing parameters, up to six solutions are predicted by the parallel flow approximation theory (four stable solutions and two unstable ones). As demonstrated in this section, the bistability phenomenon is not affected by the aspect ratio of the enclosure or by the type of boundary conditions applied on the system.

As a complementary literature review, within the context of double-diffusive convection, the existence of multiple solutions has already been reported in the past by a few authors. For instance, Kalla et al. (Reference Kalla, Mamou, Vasseur and Robillard2001) investigated the problem of thermosolutal convection in a Darcy porous layer, subject to cross-fluxes of heat and mass. It was found that, for subcritical convection conditions, up to five different solutions are possible (two of which are unstable), but these were induced by a lateral heating. Here, the effect of the horizontal fluxes of heat is similar to that of the modified Prandtl number in the present study. Due to the lack of symmetry of the thermal and solutal boundary conditions, the number of stable solutions is reduced from four to three. Similar results were reported by Mansour, Amahmid & Hasnaoui (Reference Mansour, Amahmid and Hasnaoui2008) while reconsidering the same problem upon combining the effects of thermodiffusion (Soret) with those of double diffusion. Natural convection in a fluid layer of a binary mixture induced by the combined action of a shear stress and the Soret effect was considered by Mahidjiba, Bennacer & Vasseur (Reference Mahidjiba, Bennacer and Vasseur2006). Here, again, the existence of up to five different solutions (three stable and two unstable) was demonstrated analytically and numerically. Recently, the onset of Marangoni convection in a power-law fluid under a zero-gravity environment was studied by Alloui & Vasseur (Reference Alloui and Vasseur2013). For shear-thinning fluids, subcritical convection is possible, such that, for a given range of governing parameters, five different solutions have also been reported by these authors.

6.3 Cavity aspect ratio effect

Finally, the effect of the porous layer aspect ratio on subcritical convection is illustrated in figure 16. The bifurcation diagram illustrates the evolution of $\unicode[STIX]{x1D6F9}_{0}$ as a function of the normalized Rayleigh number ( $R_{T}/R^{sup}$ ), where $R^{sup}$ represents the critical Rayleigh number for the onset of thermal convection in a confined horizontal enclosure subject to a constant heat flux. The constant $R^{sup}$ values are given in figure 16 for different values of  $A$ . The values are computed from (4.12) reported in Mamou & Vasseur (Reference Mamou and Vasseur1999). From the figure, with a chosen value of $P_{r}^{\ast }=2$ relying on previous findings, it appears clearly that the bistability convection occurs regardless of the aspect ratio value. However, it can be seen that the extent of the bistability region depends on the aspect ratio, and becomes relatively narrow when the aspect ratio is small. The region remains unaffected for large aspect ratio according to the asymptotic solution corroborated by the numerical solution of the full governing equations.

Figure 16.

Bifurcation diagram in terms of $\unicode[STIX]{x1D6F9}_{0}$ as a function of $R_{T}/R^{sup}$ for $P_{r}^{\ast }=2$ , $\unicode[STIX]{x1D711}=-0.1$ , $Le=5$ , $\unicode[STIX]{x1D700}=1$ and different values of $A$ , with $(A,R^{sup})=(1,22.94)$ , (2, 14.79), (3, 13.24), (5, 12.45) and (10, 12.11). The lower and upper curves are joined manually by dotted curves for a better display of the bifurcation region.