Nomenclature
d,	depth of the porous layer;	PrD,	Darcy-Prandtl number;
g,	gravitational acceleration;	q,	velocity vector;
K,	permeability of the porous medium;	RSi,	solute Darcy-Rayleigh number of the ith-component;
k,	unit vector in the vertical direction;	RT,	thermal Darcy-Rayleigh number;
M,	ratio of heat capacities;	t,	time;
p,	pressure;	x, y, z,	space coordinates.
Greek symbols
α,	horizontal wave number;	Λ1,	stress relaxation parameter;
αT,	thermal expansion coefficient;	Λ2,	strain retardation parameter;
αSi,	solute analogue of αT, i = 1, 2;	µ,	dynamic viscosity;
ε,	porosity;	ν,	kinematic viscosity;
κT,	thermal diffusivity;	ρ,	fluid density;
κSi,	solute diffusivity, i = 1, 2;	σ,	growth term;
λ1,	stress relaxation time;	τi,	ratio of diffusivities, i = 1, 2;
λ2,	strain retardation time;	ψ,	stream function.
Subscripts/superscripts
b,	basic state;	U,	upper boundary;
L,	lower boundary;	*,	dimensionless variable.

1 Introduction

The study of convective instabilities is one of the central problems in porous media fluid mechanics because of its natural occurrence and relevance in many science and engineering applications. It is a well-established fact that, in a horizontal porous layer, thermal instability sets in as stationary convection, and no oscillatory motion takes place^[1-2]. Nonetheless, in the case of double diffusive convection, it is intriguing to note that oscillatory convection appears as a preferred mode of instability under some conditions. Copious literatures available on these topics are documented in the books^[3-6].

Of late, researchers have been showing an increasing attention in the study of natural convection of non-Newtonian fluids saturating porous media, and the studies are very sparse. Most of the existing studies have focused on thermal convective instability in a viscoelastic fluid-saturated porous layer^[7-10] as applications of these fluids are found in geophysics, biorheology, pharmaceutical and petroleum industries, polymer processing and so on. An important result noticed that the fluid elasticity influences the preferred mode of instability from the stationary mode to the oscillatory mode. Besides, investigations have also been carried out to investigate the consequence of additional diffusing component (double diffusive convection) on thermal convective instability in a viscoelastic fluid-saturated porous layer^[11-13].

Convective instability in ternary fluids is a complex process, and the presence of thermal currents and/or concentrations leads to linear and nonlinear behaviours. In a ternary fluid, the density depends on three different stratifying agents, and one of which may be the temperature. This leads to a competition between different diffusions of stratifying agents, and as a result, oscillatory convection may occur. Triple diffusive convection is found to occur in many practical applications. For example, aqueous suspensions of deoxyribonucleic acid (DNA) may contain more than two independently diffusing components with different diffusivities. The DNA coils are bound to be extended as they are subject to the action of shear or externally applied tension, and as a result, the viscoelastic properties of DNA are to be taken into consideration in the study of DNA convection. Although the oscillatory convection in viscoelastic fluids is not observed in experimental conditions^[14], Perkins et al.^[15-16] showed that such a convection persists as the first convective instability in dilute suspensions of long DNA molecules. Later, Kolodner^[17] also confirmed this possibility. Another important situation stems from the study of medical screening tool for abnormalities in lipids, containing various components such as cholesterol (low-density lipoproteins, high-density lipoproteins) and triglycerides which possess different diffusivities and exhibit viscoelastic properties. These proteins and fats help to generate energy in the human body, and this complex process is collectively known as cellular respiration. In particular, the variations in these diffusivities help to determine approximate risks for cardiovascular and other diseases. Another area of interest is in understanding contaminant transport^[18-19], where the chemical species that form the contaminants are non-reactive. These observations have triggered a renewed interest in the study of linear and nonlinear convective instabilities in viscoelastic triple diffusive fluid systems.

The triple diffusive convection in a Newtonian fluid layer has been studied extensively^[20-27]. In particular, Pearlstein et al.^[24] discussed several applications and established some novel results which were overlooked by the previous investigators and also not found in double diffusive fluid systems. Shivakumara and Naveen-Kumar^[28] investigated linear and nonlinear multi-diffusive convection in a horizontal layer of couple stress fluids. The study of triple diffusive convection in porous media has also been analyzed. Rudraiah and Vortmeyer^[29] and Poulikakos^[30] studied the onset of triple diffusive convection in a porous medium by following the analysis of Griffiths^[20]. Whereas Tracey^[31] followed the analysis of Pearlstein et al.^[24] to study the triple diffusive convection in porous media, and Rionero^[32-33] extended the study to obtain some additional results. The linear instability of the triple diffusive Maxwell fluid-saturated porous layer was considered by Zhao et al.^[34] while a weakly nonlinear stability analysis of the problem has been investigated recently by Raghunatha et al.^[35].

The studies of convective instability in viscoelastic fluids saturating a porous medium are mainly focused either on single or double diffusive fluid systems despite the multi-diffusive convection is bound to occur in many practical situations mentioned above. As far as the Newtonian fluid systems are concerned, some fundamental differences between double and triple diffusive fluid systems in the onset and development of convective motion are unveiled. In light of these observations, it is pragmatic that to explore the qualitatively new features found by Tracey^[31] will carry over to the case of Oldroyd-B fluids saturating a porous medium as well since such a study sheds more light on the complex instability characteristics of the system. The objective of the current study is to investigate linear and weakly nonlinear triple diffusive convection in an Oldroyd-B fluid-saturated porous medium and to unveil some interesting results that the system is capable of supporting. Also, the main purpose is to understand the stability of oscillatory bifurcating solution by deriving a cubic Landau equation and also to understand the influence of various physical parameters on convective heat and mass transfer.

2 Governing equations

The physical configuration is shown in Fig. 1. We consider an incompressible Oldroyd-B type of a viscoelastic fluid-saturated horizontal layer of a Darcy porous medium of the thickness d containing three stratifying agents with different molecular diffusivities, where all the three may be solute concentrations, or one of the components may be the temperature T, and the other two are solute concentrations S_i (i=1, 2). The lower impermeable boundary is kept at the higher temperature and solute concentrations, while at the upper impermeable boundary, they are maintained at lower values. A Cartesian coordinate system (x, y, z) is chosen such that the origin is at the bottom of the porous layer. The gravity is acting in the negative vertical z-direction. A condition of local thermal equilibrium between the solid phase and the fluid phase of the porous medium is assumed.

A modified Darcy-Oldroyd model relevant to the case of a Darcy porous medium has been considered to describe the flow in a porous medium in the form^[34-38] of

(1)

where q=(u, v, w) is the velocity, p is the pressure, μ is the fluid viscosity, g is the acceleration due to gravity, K is the permeability of the medium, ρ is the fluid density, ε is the porosity, λ₁ is the relaxation time, and λ₂ is the retardation time. It may be noted that the time derivative term has been included in the above equation as it contributes for the convection to set in as oscillatory motions. The viscoelastic model considered includes the classical viscous Newtonian fluid as a special case for λ₁=λ₂=0 and the Maxwell fluid for λ₂=0. The conservation equations for the stratifying agents are

(2)

(3)

(4)

where M=(ρ₀c)_m /(ρ₀c)_f=((1-ε)(ρ₀c)_s+ε (ρ₀c)_f)/(ρ₀c)_f is the ratio of heat capacities, c is the specific heat, κ_T is the thermal diffusivity, and κ_S1 and κ_S2 are the solute analogues of κ_T. The subscripts m, f, and s refer to the porous medium, fluid, and solid, respectively. Under the Oberbeck-Boussinesq approximation, the equation of state is

(5)

where α_T is the thermal expansion coefficient, α_S1 and α_S2 are the solute analogues of α_T, and ρ₀ is the density at the reference temperature and solute concentrations.

The boundary conditions are

(6)

(7)

(8)

where S_iU < S_iL and T_U < T_L with the sign convention that ΔS_i=S_iL-S_iU>0 when a solute concentration is stabilizing and ΔT=T_L-T_U>0 when the temperature is destabilizing.

The above governing equations are made dimensionless using the following transformations:

(9)

Substituting Eq. (9) into the governing equations, we obtain

(10)

(11)

(12)

(13)

where Λ₁=λ₁κ_T/ (d²ε) is the stress relaxation time parameter, Λ₂=λ₂κ_T/ (d²ε) is the strain retardation time parameter, R_T=α_TgΔTKd /(νκ_T) is the thermal Darcy-Rayleigh number, R_Si=α_SigΔS_iKd /(νκ_T) (i=1, 2) is the solute Darcy-Rayleigh number of the ith-component, Pr_D=νd²ε²/ (κ_TK) is the Darcy-Prandtl number or the Vadasz number, τ_i=κ_Si/κ_T (i=1, 2) is the ratio of diffusivities, and A=M /ε. If all the components are solute concentrations, then the value of A turns out to be unity. The basic state is quiescent, and temperature, solute concentrations, and pressure distributions are found to be

(14)

where p₀ is the pressure at z=0, and the subscript b denotes the basic state. To investigate the instability of the basic state, we now superimpose infinitesimal perturbations on the basic state of the form,

(15)

where primes indicate the perturbed quantities. Substituting Eq. (15) into Eqs. (10)-(13), eliminating the pressure term from the momentum equation by operating curl and introducing a perturbation stream function ψ' (x, y, z) with

(16)

we obtain the governing nonlinear stability equations which can be written as

(17)

where J(·, ·) is the Jacobian function with respect to x and z, representing the nonlinear terms, and L is the linear differential operator given by

(18)

where .

The boundary conditions now become

(19)

3 Linear instability analysis

The linear instability analysis proceeds by ignoring the nonlinear terms in Eq. (17). Then, we eliminate T, S₁, and S₂ from the momentum equation to provide one equation for ψ in the form,

(20)

As usual, a normal mode solution is assumed as

(21)

where α and π are wave numbers in the x - and z -directions, respectively, and σ(=σ_r+iω) is the growth term. Substituting Eq. (21) into Eq. (20) yields the following expression:

(22)

where δ²=α²+π². To investigate the instability of the system, the real part of σ is set to zero, and let σ=iω in Eq. (22) and also clear the complex quantities from the denominator to obtain

(23)

where

(24)

Since R_T is a physical quantity, it implies either ω=0 or N=0 in Eq. (23). The condition ω=0 corresponds to the stationary convection, and the case N=0 (ω ≠ 0) corresponds to the oscillatory convection.

3.1 Stationary convection

The condition ω=0 corresponds to the stationary onset, and thus

(25)

is the thermal Darcy-Rayleigh number, where the layer is unstable with respect to the stationary onset. It is observed that Eq. (25) is independent of viscoelastic parameters, indicating that the viscoelastic effects appear only in the case of time dependent motions, and coincide with the result obtained for ordinary viscous fluids^[29]. We note that, for the stationary neutral solution, R_T^S is a single valued function of α and attains its critical value at α_c=π. Therefore, the critical thermal Darcy-Rayleigh number for the stationary onset is

(26)

In the absence of additional diffusing components, Eq. (26) gives R_Tc^S=4π², which is the known exact value.

3.2 Oscillatory convection

For the oscillatory case, ω ≠ 0 and hence N=0, which gives a dispersion relation cubic in ω²,

(27)

where

When N=0, Eq. (23) shows that the oscillatory convection occurs at R_T=R_T⁰, where

(28)

and ω² is given by Eq. (27). Equation (27) shows that, for a suitable combination of parameters A, Pr_D, Λ₁, Λ₂, τ₁, τ₂, R_S1, and R_S2, it is possible to have more than one positive values of ω² at the same wave number α, indicating the existence of the multiple oscillatory neutral solution R_T⁰ which has important implications on the nature of the instability of the system. For a single diffusive case (R_S1=R_S2=0), it is observed that the oscillatory convection occurs provided that

(29)

Thus, the necessary condition for the occurrence of oscillatory convection is

(30)

For a double diffusive case (R_S2=0, say), we derive a dispersion relation quadratic in ω²,

(31)

where

Since ω²>0 for the occurrence of oscillatory convection, a cautious glance at Eq. (31) provides the imminent conditions as

(32)

The point at which the steady and oscillatory neutral curves meet is called the bifurcation point which can be located in the (α, R_T) plane. The bifurcation points will occur on the steady neutral curve at α=α_b for which ω=0 satisfying Eq. (27). Thus, a₄ (α_b)=0, or equivalently

(33)

where φ_i=Pr_D(Aτ_i-1) (i=1, 2) and ϖ=(Λ₂-Λ₁) Pr_D+1. The minimum value of R_T⁰ with respect to α, denoted by R_Tc⁰, is calculated using the procedure explained by Raghunatha et al.^[35].

4 Oscillatory nonlinear stability analysis

The linear instability theory discussed in the previous section gives us the condition for instability. Here, we use the perturbation method to perform the nonlinear stability analysis and derive the cubic Landau equation. Such a study helps to analyze the convective rate of heat and mass transfer. Since the viscoelastic effects are time dependent, the study has been restricted to the periodic convective solution that bifurcates at R_T=R_Tc⁰. Accordingly, a small parameter χ that indicates the deviation from the critical point is introduced. Consequently, the stream function, temperature, and solute concentrations are expanded in terms of χ as^[11]

(34)

The fast time scale t and the slow time scale s are introduced in the form of s=χ² t, and hence the operator is replaced by the operator . Substituting Eq. (34) into Eq. (17), we get a set of equations at each order of χ. At the leading order in χ, the coupled differential equation is linear and homogeneous,

(35)

which corresponds to the linear instability problem for the oscillatory convection discussed in the previous section. The eigenvalue is R_T=R_Tc⁰, and the eigenfunctions are of the form

(36)

(37)

(38)

(39)

where the overbar denotes the complex conjugate. The amplitudes A₁-D₁ and A₁-D₁ are functions of the slow time s, whereas ω and α are taken to be the critical values associated with R_T=R_Tc⁰. The undetermined amplitudes are related by

(40)

The inhomogeneous terms appearing in the equations at the order χ² are evaluated using the first order solutions,

(41)

(42)

(43)

From the above relations, we can deduce that the corrections to the stream function, temperature, and solute concentrations fields have terms of the frequency 2ω. Thus, the second order stream function, temperature, and solute concentrations fields can be expressed as follows:

(44)

(45)

(46)

(47)

where (ψ₂₀, ψ₂₂), (T₂₀, T₂₂), and (S₁₂₀, S₁₂₂, S₂₂₀, S₂₂₂) are, respectively, the stream functions, the temperatures, and solute concentrations fields which have terms of the frequency 2ω and independent of the fast time scale t.

At the second order in χ, the equation is

(48)

The above equation is now solved, and the solutions are found to be

(49)

At the third order in χ, the equation now becomes

(50)

where

The inhomogeneous terms appearing in Eq. (50) are evaluated using the first and second order solutions, and the explicit expressions are given in Appendix A. The third order equations have the solutions

(51)

(52)

(53)

(54)

Applying the solvability condition to Eq. (50) yields the following cubic Landau equation:

(55)

where

(56)

(57)

with

Writing B₁ in the phase-amplitude form,

(58)

and substituting this into Eq. (55), we get the following expressions for the amplitude |B₁| :

(59)

(60)

where , γ^-1η=l_r+il_i, and p_h(·) represents the phase shift. The magnitude and direction of the periodic convective solution and the frequency shift are obtained in Eq. (55). Our curiosity is in the direction of the bifurcation, which is easily shown to be dependent on the sign of the quantity,

(61)

If Q is positive, the bifurcation is supercritical and stable, and it is subcritical and unstable if Q is negative. For this case, the time and area-averaged thermal Nusselt number (Nu_T) and Solute Nusselt numbers (Nu_S1, Nu_S2) are determined, and they are given by

(62)

(63)

5 Results and discussion

Both linear instability and weakly nonlinear stability of an initially quiescent triply diffusive Oldroyd-B viscoelastic fluid-saturated porous layer are investigated. There are seven nondimensional parameters Pr_D, Λ₁, Λ₂, τ₁, τ₂, R_S1, and R_S2, which influence the stability characteristics of the system, and the thermal Darcy-Rayleigh number is taken as a free parameter. The results are presented for a wide range of the known physical parameters. The value of Pr_D is allowed to vary from 10 to 200 (see Ref. [39]). The viscoelastic parameters Λ₁ and Λ₂ are chosen to be either greater or less than unity ^[40]. In particular, the results are also presented for the case of low-density lipoproteins, high-density lipoproteins and triglycerides system for which the diffusion coefficients are, respectively^[41], κ_T=0.35×10^-3 cm²/s, κ_S1=0.46×10^-4 cm²/s, and κ_S2=0.53×10^-4 cm²/s. The presence of a third diffusing component supports some unusual dynamical behaviours which are not observed in the study of single/double diffusive viscoelastic fluid systems. Some novel results have been delineated by comparing the linear instability characteristics of Oldroyd-B, Maxwell, and Newtonian fluids as particular cases.

5.1 Neutral stability curves

The neutral stability curves in (α, R_T⁰) and (α, ω²) planes for different values of the strain retardation parameter Λ₂, the Darcy-Prandtl number Pr_D, and the stress relaxation parameter Λ₁ are displayed in Figs. 2-5. It is observed that there exists only one positive value of ω² for the parametric values chosen in these figures. The oscillatory neutral stability curves in the (α, R_T⁰) plane show an upward concave shape, and the region below each such curve confines to the region of stability, while the region above it corresponds to instability. It is seen that the oscillatory convection occurs even if the ratio of diffusivities is greater than unity, and a reverse result is noticed compared with double diffusive fluid systems (see Fig. 3(a)). Besides, increasing Λ₂ is to increase the value of R_T⁰, indicating that its effect is to enhance the region of stability and also to suppress the frequency of oscillations (see Figs. 2(b) and 3(b)). Moreover, increasing Pr_D (see Fig. 4(a)) and Λ₁ (see Fig. 5(a)) is to decrease the region of stability with the tendency of increasing the frequency of oscillations (see Figs. 4(b) and 5(b)). The results shown in Figs. 6(a) and 6(b) for the case of low-density lipoproteins, high-density lipoproteins and triglycerides system exhibit a similar kind of behaviour as seen in Figs. 2(a) and 2(b). The stabilizing/destabilizing influence of an additional diffusing component on neutral stability curves is displayed in Fig. 7(a) by taking positive/negative values of R_S1, respectively, for Pr_D=50, A=2, Λ₁=0.3, Λ₂=0.2, τ₁=0.3, τ₂=0.4, and R_S2=100. It is observed that the region of stability increases when the diffusing component is more stabilizing, and also its effect is to suppress the frequency of oscillations (see Fig. 7(b)).

5.2 Disconnected oscillatory neutral curves

The existence of disconnected oscillatory neutral curves under certain conditions is an integral feature of triple diffusive fluid systems which has an important consequence on the linear instability of the system. Such a situation can be obtained by locating the bifurcation points. For this, Eq. (33) has to be solved for any chosen parametric values, which shows that there is a possibility of getting either zero or one/two/three bifurcation points.

The viscoelastic parameters Λ₁ and Λ₂ can be greater than unity for many polymeric fluids, and it is interesting to discern the evolution of neutral stability curves for this case as well. Figures 8(a)-8(f) exhibit evolutions of neutral stability curves for various positive values of R_S2 ranging from 66.71 to 68 for Pr_D=50, A=1.7, Λ₁ =1.32, Λ₂=1.2 (Λ₁, Λ₂>1), τ₁=0.191, τ₂= 0.264, and R_S1=-30. Figure 8(a) shows that the oscillatory neutral curve is joined to the stationary neutral curve at two bifurcation points. As the value of R_S2 goes on decreasing, the bifurcation points move closer together and eventually detach from the stationary neutral curve (see Figs. 8(b) and 8(c)). Figures 8(b) and 8(c) also show that the left side upper maximum of the closed oscillatory neutral curve lies above the minimum of the steady neutral curve, suggesting that the instability of the system can be conveniently analyzed by a single critical value of R_T. Figure 8(d) shows the complete detachment of the closed convex oscillatory neutral curve from the stationary curve with further decrease in the value of R_S2, and when R_S2=66.75, the oscillatory curve lies well below the minimum of the stationary neutral curve as shown in Fig. 8(e). Figure 8(e) clearly demonstrates that three critical values of R_T are needed to identify the instability criteria instead of the usual single value. It is seen that the layer is linearly stable if R_T < R_T1 and R_T2 < R_T < R_T3, while for R_T1 < R_T < R_T2 and R_T>R_T3, the layer is unstable. The closed oscillatory neutral curve finally collapses to a point and disappears leaving only the stationary neutral curve when R_S2=66.71. It is important to note here that the closed disconnected oscillatory neutral curve is not perfectly heart-shaped, a result that is in contradiction to the Newtonian fluid case^[31] in which the disconnected oscillatory neutral curve has two extrema at the same thermal Darcy-Rayleigh number (i.e., quasiperiodic bifurcation). Thus, we note that this degeneracy is broken due to the presence of viscoelastic effects which contribute in an opposite way on the onset of oscillatory convection.

The critical thermal Darcy-Rayleigh numbers depicting the multivalued region are exhibited in Fig. 9 for the case considered in Fig. 8. To the right of the point A of the infinite slope (R_S2>67.075), oscillatory instability first occurs at a lower value of R_Tc than the stationary instability, and there is a single value of R_Tc. To the left of the point B (R_S2 < 66.724 2), instability occurs as the stationary convection, and oscillatory instability does not occur, and again there is one value of R_Tc. However, for the finite range of values of R_S2 lying between the vertical lines passing through the points A and B, three critical values of the thermal Darcy-Rayleigh number are needed to delineate the linear instability criteria.

Figures 10(a)-10(f) represent evolutions of neutral curves for different values of R_S2 when Pr_D=30, A=1.85, Λ₁=0.273, Λ₂ =0.173 (Λ₁, Λ₂ < 1), τ₁=0.19, τ₂=0.273, and R_S1=-22.50. The behaviour of neutral curves is similar to that observed in the previous case, but the range of R_S2 differs in which the disconnected oscillatory neutral curves occur. In this case, the closed convex oscillatory neutral curve goes on dwindling and eventually disappears leaving only the stationary neutral curve at R_S2=55.3 (see Fig. 10(f)). Moreover, in this case, the disconnected oscillatory neutral curve is found to be not heart-shaped. The stability boundaries correspond to this case, as shown in Fig. 11.

It is intriguing to note that the above observed phenomena also carry over to the case of low-density lipoproteins, high-density lipoproteins and triglycerides system as well, and the same phenomenon is evident from Figs. 12(a)-12(f). The results presented here are for Pr_D=10, A=1, Λ₁=0.2, Λ₂=0.1 (Λ₁, Λ₂ < 1), τ₁=0.13, τ₂=0.15, and R_S1=-10. These figures indicate that the results are similar to those for the other transport property ratios. From Fig. 12(e), we observe that, the closed convex oscillatory neutral curve has become topologically detached from the stationary neutral curve, demonstrating the requirement of three values of critical thermal Darcy-Rayleigh number to specify the linear instability criteria under certain conditions. Finally, it is seen that the closed convex oscillatory neutral curve collapses to a point and vanishes, leaving only the stationary neutral curve. The corresponding stability boundaries are presented in Fig. 13.

Figures 14 and 15 accomplish the sensitivity of viscoelastic parameters Λ₁ and Λ₂ on the instability characteristics of the system, respectively. Figures 14(a) and 15(a) suggest three values of critical thermal Darcy-Rayleigh number to identify the instability of the system. However, a minor variation in the value of either Λ₁ (see Fig. 14(b)) or Λ₂ (see Fig. 15(b)) completely changes the mode of instability itself. That is to say, the mode of instability switches over from oscillatory to stationary, indicating the adequacy of a single critical thermal Darcy-Rayleigh number to specify the linear instability criteria.

The similarities and differences of the predictions of the two viscoelastic models as well as the Newtonian fluid are enunciated in Figs. 16(a)-16(c). Figures 16(a), 16(b), and 16(c), respectively, display the nature of neutral stability curves for Oldroyd-B (Λ₁=0.15, Λ₂=0.1), Maxwell (Λ₁=0.15, Λ₂=0) and Newtonian (Λ₁=Λ₂=0) fluids with Pr_D=30, A=2, τ₁=0.201, τ₂=0.249, R_S1=-23, and R_S2=48.5. In the case of an Oldroyd-B fluid, three critical thermal Darcy-Rayleigh numbers are required to specify the linear instability criteria, while in the case of remaining two types of fluids, a single critical thermal Darcy-Rayleigh number is enough to describe the linear instability. Moreover, the onset of instability ceases to be oscillatory in the case of Newtonian fluids, and it is seen that the instability sets in as only stationary convection.

5.3 Critical condition

Analyzing the behaviour of critical thermal Darcy-Rayleigh number R_Tc and the corresponding critical frequency of oscillations ω_c for various physical parameters is important as it provides the crucial information on the linear instability of the system. The variations of R_Tc and ω_c as a function of Λ₂ are shown in Figs. 17(a) and 17(b), respectively, for different values of Λ₁ when A=1.78, τ₁=0.208, τ₂=0.244, R_S1=-20, and R_S2=30. Figure 17(a) clearly exposes the dual effect of elasticity parameters on the nature of onset of instability. It is seen that, for a fixed value of Λ₁, there exists a threshold value of Λ₂ below which the instability first sets in as oscillatory convection and beyond which stationary convection is preferred. Moreover, the threshold value of Λ₂ increases with increasing Λ₁. The onset of oscillatory convection gets delayed with increasing Λ₂ while the trend gets reversed with increasing Λ₁. This is because the increase in Λ₂ amounts to the increase in the time taken by the fluid element to respond to the applied stress. That is, during the growth of the retardation parameter, the effect of friction will be higher, and as a result higher values of R_Tc are needed to instil instability on the system with increasing Λ₂. Whereas the increase in Λ₁ amounts to allowing the applied stress to act for a longer time on the fluid element which results in lesser elastic memory of the fluid. In fact, the increase of relaxation ceases the stickiness of the fluid, and hence the effect of friction will be lesser so that the convection sets in at lower values of R_Tc with increasing Λ₁. The opposite behaviour of elasticity parameters on the onset of oscillatory convection is also seen through the variation in ω_c shown in Fig. 17(b). From this figure, it is observed that ω_c increases with increasing Λ₁ due to an increase in the elasticity of the fluid, and an opposite behaviour could be seen with increasing Λ₂. Thus, it is evident that the increase in the overstability vibration is to advance the onset of oscillatory convection. In Fig. 18, we display the effect of Darcy-Prandtl number Pr_D on the critical thermal Darcy-Rayleigh number R_Tc. We observe that the effect of increasing Pr_D is to increase the range of threshold value of Λ₂ at which the instability switches over from oscillatory to stationary and also to hasten the onset of oscillatory convection.

It is possible to estimate R_Tc for both stationary and oscillatory convections which coincide for some parametric values, and as a result, a codimension-two bifurcation occurs. This behaviour is illustrated in Fig. 19 on the (Λ₂/Λ₁, Λ₁) plane for different known values of Pr_D obtained from Eqs. (26) and (28). In the figure, each curve locates the boundary separating stationary and oscillatory solution regions. For fixed values of Λ₂/Λ₁ and Pr_D, there exists a value of Λ₁=Λ₁^*, where R_Tc^S=R_Tc⁰. The system is unstable under oscillatory convection in the region above each curve, and below the curve, the system is unstable under stationary convection. For a fixed value of Λ₂/Λ₁, the value of Λ₁ at which a codimension-two bifurcation occurs decreases with increasing Pr_D.

5.4 Oscillatory nonlinear stability

The stability of oscillatory bifurcating solutions is analyzed by deriving a cubic Landau equation. The nature of oscillatory bifurcation can be determined from the sign of Q, and the point at which Q changes the sign is termed as the point. Figures 20-22 represent the computed values of Q as a function of R_S2 for different values of physical parameters, namely, Pr_D, R_S1, Λ₁, Λ₂, τ₁, and τ₂. These figures demonstrate the possibility of subcritical oscillatory bifurcation for a range of parametric values, and indicate the occurrence of instability before the linear threshold is reached. This is expected, because the linear instability analysis only provides a sufficient condition for instability. Besides, the tricritical point is shifted towards higher values of R_S2 with increasing Λ₁, R_S1, and τ₂, while an opposite trend could be seen with increasing Pr_D, Λ₂, and τ₁.

5.5 Heat and mass transport

Heat and mass transfer are measured through thermal and solute Nusselt numbers. The thermal Nusselt number (Nu_T) and solute Nusselt numbers (Nu_S1, Nu_S2) are displayed as a function of R_T for different values of Pr_D, Λ₁, and Λ₂ in Figs. 23(a), 23(b), and 23(c), respectively. The increase in the value of Λ₁ is to increase the heat and mass transfer characteristics of oscillatory convection, whereas the increase in Λ₂ and Pr_D exhibits an opposite kind of behaviour. This may be attributed to the fact that the increase in Λ₁ is to increase the overstability vibration (see Figs. 5(b) and 17(b)), and consequently, its effect is to enhance the heat and mass transfer, while the increase in Λ₂ is to suppress the heat and mass transfer because their effect is to decrease the overstability vibration (see Figs. 3(b) and 17(b)).

6 Conclusions

The triple diffusive convection in an Oldroyd-B fluid-saturated porous medium is investigated. By performing the linear instability analysis, the condition for the onset of stationary and oscillatory convection is established. It is found that the oscillatory convection occurs even if the ratio of diffusivities is greater than unity. Under certain conditions, disconnected closed convex oscillatory neutral curves are found to occur, indicating the requirement of three critical thermal Darcy-Rayleigh numbers to specify the linear instability criteria instead of the usual single value. This result is in contradiction to the observed phenomenon for Newtonian fluids wherein the disconnected oscillatory neutral curve is found to be perfectly heart-shaped (i.e., quasiperiodic bifurcation). The effect of increasing retardation parameter is to delay the onset of oscillatory convection, while the trend gets reversed with the increasing relaxation parameter. The value of retardation parameter beyond which the instability sets in as oscillatory convection increases with the increasing relaxation parameter as well as the Darcy-Prandtl number. The instability characteristics of the system analyzed under the same parametric values for Oldroyd-B, Maxwell, and Newtonian fluids are found to be qualitatively different. A perturbation method is used to perform a weakly nonlinear stability analysis, and subcritical bifurcation is found to be possible depending on the choice of physical parameters. Also, the increase in the value of the relaxation parameter is to enhance the time and area-averaged heat and mass transfer, while the increase in the retardation parameter and the Darcy-Prandtl number exhibits an opposite kind of performance.

Appendix A  

The inhomogeneous terms appearing in the third order equations are calculated using the first and second order solutions,

Acknowledgements The authors thank reviewers for their constructive remarks and useful suggestions, which improve the work significantly.