1. Introduction
The pioneering paper by Gill (Reference Gill1969) has provided an interesting analytical proof that convection due to thermal buoyancy in a vertical porous layer bounded by isothermal impermeable walls is always stable for all infinitesimal perturbations. The obtained results offered significant justification for the use of insulating porous materials in buildings, obviating the necessity for air gaps. The work by Gill (Reference Gill1969) was further extended by other authors including additional effects such as the Prandtl–Darcy number (Rees Reference Rees1988), lack of local thermal equilibrium between the fluid and the solid phases (Rees Reference Rees2011) and the nonlinearity of the perturbation dynamics (Straughan Reference Straughan1988; Scott & Straughan Reference Scott and Straughan2013). These investigations collectively furnished compelling evidence that the fluid flow stability established by Gill remains valid. On the other hand, the possibility of Gill's problem becoming unstable was established by taking into account temperature-dependent viscosity (Kwok & Chen Reference Kwok and Chen1987; Shankar, Nagamani & Shivakumara Reference Shankar, Nagamani and Shivakumara2023), vertical boundaries to be permeable (Barletta Reference Barletta2015), boundary and inertia effects (Shankar, Kumar & Shivakumara Reference Shankar, Kumar and Shivakumara2017), viscoelasticity of the fluid (Shankar & Shivakumara Reference Shankar and Shivakumara2017; Lazzari, Celli & Barletta Reference Lazzari, Celli and Barletta2021), maximum density property with time-dependent velocity term in Darcy's law (Naveen, Shankar & Shivakumara Reference Naveen, Shankar and Shivakumara2020), horizontal heterogeneity in permeability (Shankar & Shivakumara Reference Shankar and Shivakumara2022), sandwiched vertical porous slab (Barletta et al. Reference Barletta, Celli, Lazzari and Brandão2022) and the Prandtl–Darcy number with internal heating (Nagamani, Shankar & Shivakumara Reference Nagamani, Shankar and Shivakumara2023).
Numerous fluid dynamical systems manifest in nature and various scenarios wherein buoyancy forces instigate the movement of a fluid laden with dissolved minerals through a porous medium. Notably, double-diffusive convection encompasses two interacting species that influence the distribution of fluid density. For instance, this phenomenon might involve two solute concentrations with distinct molecular diffusivities or combinations of heat and solute concentration. The exploration of double-diffusive convection in a fluid-saturated porous medium has remained a vibrant and enduring research domain for numerous decades, owing to its significance and broad-ranging applications across disciplines including oceanography, migration of pollutants through saturated soil, geophysics, astrophysics and engineering (Griffiths Reference Griffiths1981; Pritchard & Richardson Reference Pritchard and Richardson2007; Straughan Reference Straughan2015). Comprehensive assessments of the theoretical and empirical analyses pertaining to double-diffusive convection in a porous medium can be found in Cheng (Reference Cheng1979) and Nield & Bejan (Reference Nield and Bejan2017). However, it is noteworthy that investigation into the stability of double-diffusive convection in a vertical porous layer has been relatively limited within the existing literature.
The stability of stationary convective flow in a vertical porous layer saturated with a binary mixture by imposing different temperatures and zero mass flux at the boundaries was considered by Gershuni, Zhukhovitskii & Lyubimov (Reference Gershuni, Zhukhovitskii and Lyubimov1976). It has been shown that in the presence of sufficient longitudinal stratification, the flow becomes unstable against thermal–concentration perturbations. Khan & Zebib (Reference Khan and Zebib1981) furthered this investigation by extending it across all ranges of the salinity Darcy–Rayleigh number. Shankar, Naveen & Shivakumara (Reference Shankar, Naveen and Shivakumara2022) delved into the effects of a second diffusing component on the stability of thermal convection in a vertical porous layer. This was done under conditions involving constant but distinct temperatures and solute concentrations at impermeable vertical boundaries. Their findings revealed the formation of closed loops in the neutral stability curves and the enlargement of the instability region with increasing solute/thermal Darcy–Rayleigh number as well as the Lewis number. Overall, the study demonstrated that double-diffusive fluid systems exhibit a significantly more intricate dynamic behaviour in comparison with the conventional single-component system (Gill Reference Gill1969).
The possibility of fluid dynamical systems in porous media, where fluid density is influenced by three or more stratifying agencies with different molecular diffusivities, cannot be ruled out and, as a result, one can expect multicomponent convection. Degens et al. (Reference Degens, von Herzen, Wong, Deuser and Jannasch1973) highlighted the occurrence of double-diffusive convection in geothermal systems, where salinity of comparable concentrations of many salts can trigger convection encompassing multi-diffusing components. Another pertinent example is in contaminant transport scenarios (Celia, Kindred & Herrera Reference Celia, Kindred and Herrera1989; Chen et al. Reference Chen, Cunningham, Ewing, Peralta and Visser1994), where various non-reactive chemical species contribute to contaminant formation. Also, the effects of dyes or small temperature gradients on the accuracy of laboratory experiments on double-diffusive convection in porous media may introduce a third property that affects the density of the fluid. Furthermore, investigating convective instability in mushy layers formed during the solidification of multicomponent alloys holds significance (Anderson & Worster Reference Anderson and Worster1996; Guba & Worster Reference Guba and Worster2006). A mushy layer, akin to a porous medium (Phillips Reference Phillips1991), facilitates interstitial melt flow, underscoring the importance of studying multicomponent convection in such contexts. These intricate fluid dynamical systems are also observed in geothermally heated lakes, magmas and their laboratory models. To this end, triple-diffusive convection in a horizontal porous layer has been undertaken by a multitude of researchers (Rudraiah & Vortmeyer Reference Rudraiah and Vortmeyer1982; Poulikakos Reference Poulikakos1985; Tracey Reference Tracey1996; Rionero Reference Rionero2012, Reference Rionero2013; Raghunatha, Shivakumara & Shankar Reference Raghunatha, Shivakumara and Shankar2018; Shivakumara & Raghunatha Reference Shivakumara and Raghunatha2022 and references therein). Several remarkable departures from those of double-diffusive fluid systems have been found when a third diffusing component is present, particularly for the onset of convection.
In this context, it is pragmatic to contribute further novel findings to the present knowledge on the stability of double-diffusive natural convection in a vertical porous layer by considering a third diffusing component. To the best of our knowledge, such study has not received any attention in the literature. The aim of this paper is to investigate the linear stability of buoyancy-driven convection in a triple-diffusive fluid-saturated vertical porous layer by assuming small-amplitude perturbations of the basic flow. The resulting eigenvalue problem governing the stability of the fluid flow is numerically solved. The topology of neutral curves is outlined systematically and the stability boundary is noted to be significantly more intricate depending on the transport property ratios. It is established that the triple-diffusive fluid systems support several remarkable departures from what occurs in single- and double-diffusive systems. It is observed that more than two critical Darcy–Rayleigh numbers correspond to the linear onset of instability under certain conditions.
2. Mathematical formulation
The stability of buoyancy-driven convection in a triple-diffusive fluid saturating a vertical layer of a Darcy porous medium of width 
$2h$ is considered as shown in figure 1. The concentration 
${C_i}\textrm{ }(i = 1,2,3)$ of each individual diffusing component (one of the diffusing mechanisms may be heat) along the right- and left-hand vertical impermeable boundaries is held at different constant value 
${C_{ir}}$ and 
${C_{il}}$, respectively. The x axis is horizontal and perpendicular to the vertical planes, the z axis is pointing vertically upward and the y axis is horizontal and parallel to the planes. The gravitational acceleration 
$\boldsymbol{g}$, with modulus g, acts in the downward vertical direction. The porous medium exhibits isotropic and homogeneous characteristics, with negligible influence of viscous dissipation and the absence of volumetric heating. The local thermal equilibrium between the solid and the fluid phases is assumed and the Oberbeck–Boussinesq approximation holds. Besides, the off-diagonal (Soret, Dufour and cross-diffusion) contributions to the fluxes of the stratifying agencies are neglected. The density 
$\rho$ of an incompressible fluid is assumed to vary linearly in the form
\begin{equation}\rho = {\rho _0}\left[ {1 - \sum\limits_{i = 1}^3 {{\alpha_i}({C_i} - {C_{i0}})} } \right],\end{equation}
where 
${\rho _0}$ is the reference density at reference solute concentrations 
${C_{i0}} = ({C_{il}} + {C_{ir}})/2$ and the coefficients 
${\alpha _i}$ are constants. The governing equations under the above assumptions are
\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{q} = 0,\end{gather}
\begin{gather}\boldsymbol{\nabla }p =- {\rho _0}\sum\limits_{i = 1}^3 {{\alpha _i}({C_i} - {C_{i0}})\boldsymbol{g}} - \frac{\mu }{K}\boldsymbol{q},\end{gather}
\begin{gather}{\varepsilon _i}\frac{{\partial {C_i}}}{{\partial t}} + (\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{\nabla }){C_i} = {\kappa _i}{\nabla ^2}{C_i}\quad (i = 1,2,3),\end{gather}
where 
$\boldsymbol{q} = (u,v,w)$ is the seepage velocity, p is the pressure, 
$\mu$ is the fluid viscosity, K is the permeability, 
${\varepsilon _i} = \varepsilon = \textrm{const}\textrm{. }(i = 1,2,3)$ is the porosity, t is the time and 
${\kappa _i}\textrm{ }(i = 1,2,3)$ are the solutal diffusivities. If one of the diffusing components is heat, say the component 
$i = 1$, then 
${\varepsilon _1}$ in (2.4) stands for the ratio of heat capacities defined as 
${\varepsilon _1} = [\varepsilon {(\rho c)_f} + (1 - \varepsilon ){(\rho c)_s}]/{(\rho c)_f}$, where 
${(\rho c)_f}$ and 
${(\rho c)_s}$ are the heat capacities of the fluid and the solid matrix, respectively, and its value is assumed to be 1. That is, the effect of the heat capacity of the solid matrix is neglected and this simplification is valid when the solid matrix is composed of low-heat-capacity materials such as metals.
A sketch of the physical model.
The boundary conditions are
\begin{equation}\boldsymbol{q}\boldsymbol{\cdot }\hat{i} = 0\ \textrm{at}\ x ={\pm} h,\quad {C_i} = {C_{il}}\ \textrm{at}\ x =- h,\quad {C_i} = {C_{ir}}\ \textrm{at}\ x = h.\end{equation}The quantities are made dimensionless as follows:
\begin{gather}\begin{array}{@{}cc@{}}{\boldsymbol{q}^ \ast }
= \displaystyle\frac{\boldsymbol{q}}{{({\kappa _1}/h)}},\quad {t^ \ast }
= \frac{t}{{(\varepsilon {h^2}/{\kappa _1})}},\quad {\nabla
^ \ast } = h\boldsymbol{\nabla },\quad {p^ \ast } =
\frac{p}{{({\kappa _1}\mu /K)}},\notag\\
C_i^ \ast = \displaystyle\frac{{{C_i} - {C_{i0}}}}{{({C_{ir}} - {C_{il}})}}\quad (i
= 1,2,3).\end{array}\end{gather}Using (2.6), (2.2)–(2.5) become (after ignoring the asterisks for simplicity)
\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{q} = 0,\end{gather}
\begin{gather}\boldsymbol{\nabla }p = \sum\limits_{i = 1}^3 {{R_i}{C_i}\hat{k} - \boldsymbol{q}} ,\end{gather}
\begin{gather}\frac{{\partial {C_i}}}{{\partial t}} + (\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{\nabla }){C_i} = \frac{1}{{L{e_i}}}{\nabla ^2}{C_i}\quad (i = 1,2,3),\end{gather}
\begin{gather}\boldsymbol{q}\boldsymbol{\cdot }\hat{i} = 0,\quad {C_i} ={\pm} {1/2}\ \textrm{at}\ x ={\pm} 1\textrm{ }(i = 1,2,3).\end{gather}
Here, 
${R_i} = {\alpha _i}g({C_{ir}} - {C_{il}})hK/{\kappa _1}\nu$ and 
$L{e_i} = {\kappa _1}/{\kappa _i}$ are the Darcy–Rayleigh numbers and the Lewis numbers, respectively, for 
$i = 1,2,3$, while 
$\nu = \mu /{\rho _0}$ is the kinematic viscosity of the fluid and 
$\hat{k}$ is unit vector along the z axis. Given the evident applicability of Squire's theorem (Barletta Reference Barletta2015; Shankar, Shivakumara & Naveen Reference Shankar, Shivakumara and Naveen2021), it has prompted the introduction of the stream function 
$\psi (x,z,t)$ in the form 
$u =- \partial \psi /\partial z$; 
$w = \partial \psi /\partial x$ and consequently the pressure is eliminated by operating curl on the momentum equation, and thus we obtain
\begin{gather}\left( {\frac{{{\partial^2}\psi }}{{\partial {x^2}}} + \frac{{{\partial^2}\psi }}{{\partial {z^2}}}} \right) = \sum\limits_{i = 1}^3 {{R_i}} \frac{{\partial {C_i}}}{{\partial x}},\end{gather}
\begin{gather}\frac{{\partial {C_i}}}{{\partial t}} - \frac{{\partial \psi }}{{\partial z}}\frac{{\partial {C_i}}}{{\partial x}} + \frac{{\partial \psi }}{{\partial x}}\frac{{\partial {C_i}}}{{\partial z}} = \frac{1}{{L{e_i}}}{\nabla ^2}{C_i}\quad (i = 1,2,3).\end{gather}The boundary conditions now become
\begin{equation}\psi = 0,\quad {C_i} ={\pm} 1/2\ \textrm{at}\ x ={\pm} 1\quad (i = 1,2,3).\end{equation}2.1. Basic state
The base flow is assumed to be fully developed, unidirectional, time-independent and laminar. These assumptions reduce (2.11) and (2.12) to ordinary differential equations of the form
\begin{gather}\frac{{{\textrm{d}^2}{\psi _b}}}{{\textrm{d}\kern0.06em {x^2}}} = \sum\limits_{i = 1}^3 {{R_i}\frac{{\textrm{d}{C_{ib}}}}{{\textrm{d}\kern0.06em x}}} ,\end{gather}
\begin{gather}\frac{{{\textrm{d}^2}{C_{ib}}}}{{\textrm{d}\kern0.06em {x^2}}} = 0\quad (i = 1,2,3),\end{gather}where the suffix b serves to denote ‘basic flow’. The associated boundary conditions are
\begin{equation}{\psi _b} = 0,\quad {C_{ib}} ={\pm} 1/2\ \textrm{at}\ x ={\pm} 1\quad (i = 1,2,3).\end{equation}The solution to the basic flow is found to be
\begin{gather}{\psi _b} = \frac{{({R_1} + {R_2} + {R_3})}}{4}({x^2} - 1),\end{gather}
\begin{gather}{C_{ib}} = \frac{x}{2}\quad (i = 1,2,3).\end{gather}
It may be noted that the nature of the basic velocity profile is not altered due to the consideration of additional solute concentration fields but the additional stratifying agencies change the magnitude of 
${\psi _b}$.
2.2. Disturbance equations
In order to test the stability of the basic flow, the standard procedure is perturbing the system. Accordingly, disturbances of small amplitude are superimposed in the form
\begin{equation}\psi = {\psi _b} + \psi ^{\prime},\quad {C_i} = {C_{ib}} + {C^{\prime}_i}\quad (i = 1,2,3),\end{equation}where primes serve to denote the perturbation field. Substituting equation (2.19) into (2.11)–(2.13) and linearizing we get (after ignoring the primes for simplicity)
\begin{gather}\left( {\frac{{{\partial^2}\psi }}{{\partial {x^2}}} + \frac{{{\partial^2}\psi }}{{\partial {z^2}}}} \right) = \sum\limits_{i = 1}^3 {{R_i}} \frac{{\partial {C_i}}}{{\partial x}},\end{gather}
\begin{gather}\frac{{\partial {C_i}}}{{\partial t}} - D{C_{ib}}\frac{{\partial \psi }}{{\partial z}} + D{\psi _b}\frac{{\partial {C_i}}}{{\partial z}} = \frac{1}{{L{e_i}}}{\nabla ^2}{C_i}\quad (i = 1,2,3),\end{gather}
\begin{gather}\psi = 0 = {C_i}\ \textrm{at}\ x ={\pm} 1\quad (i = 1,2,3),\end{gather}
where the operator 
$D = \textrm{d}/\textrm{d}\kern0.06em x$.
2.3. Normal modes
We now focus on the dynamics of normal modes expressed as
\begin{equation}(\psi ,{C_i})(x,z,t) = (\varPsi ,{\varPhi _i})(x)\,{\textrm{e}^{\textrm{i}a(z - ct)}},\end{equation}
with a real wavenumber a and temporal complex eigenvalue 
$c\textrm{ }( = {c_r} + \textrm{i}{c_i})$. We now substitute equation (2.23) into (2.20)–(2.22), which yields the eigenvalue problem
\begin{gather}({D^2} - {a^2})\varPsi = \sum\limits_{i = 1}^3 {{R_i}D{\varPhi _i}} ,\end{gather}
\begin{gather}- \textrm{i}aD{C_{ib}}\varPsi + \textrm{i}aD{\psi _b}{\varPhi _i} - \frac{1}{{L{e_i}}}({D^2} - {a^2}){\varPhi _i} = \textrm{i}ac{\varPhi _i}\quad (i = 1,2,3),\end{gather}
\begin{gather}\varPsi = {\varPhi _i} = 0\textrm{ }(i = 1,2,3)\ \textrm{at}\ x ={\pm} 1.\end{gather}
The imaginary part of the eigenvalue, i.e. the growth rate 
${c_i}$, is an important parameter as it allows one to detect the linear stability 
$({c_i} \le 0)$ or instability 
$({c_i} > 0)$ of the basic solution. We could not find a rigorous proof that the solution of (2.24)–(2.26) yields 
${c_i} \le 0$, as proved by Gill (Reference Gill1969) for a single-diffusing-component system. We may expect that the instability initiates when the parameter 
${R_1}$ exceeds a threshold value by considering the other governing parameters as given. Thus, the above eigenvalue problem allows one to evaluate numerically the threshold value of 
${R_1}$ for the onset of instability.
3. Numerical procedure
To solve the stability eigenvalue problem governed by (2.24)–(2.26), the spectral Chebyshev collocation technique has been employed. Accordingly, the governing equations are discretized along the x direction at the Gauss–Lobatto collocation points of 
$M\textrm{th}$ order in the domain 
$[ - 1,1]$. The Gauss–Lobatto points are given as
\begin{equation}{x_j} = \cos \left( {\frac{{{\rm \pi}j}}{M}} \right),\quad j = 0,1,2, \ldots ,M,\end{equation}
where M represents the order of the base polynomial. These 
$M + 1$ grid points correspond to the extreme points of the first kind Chebyshev polynomial of degree 
$M$:
\begin{equation}{\tau _M}(x) = \cos (M{\cos ^{ - 1}}x).\end{equation}
Using a truncated Chebyshev series, we construct an interpolating polynomial in terms of unknown field variables at the collocation points. The unknown function 
$f(x)$ and its derivatives are then approximated in accordance with the values of this function at the collocation points. The 
$M\textrm{th}$-order polynomial representation of function 
$f(x)$ is presented as follows:
\begin{equation}f(x) = \sum\limits_{j = 0}^M {{h_j}(x)f({x_j})} ,\end{equation}
where the Lagrange interpolant 
${h_j}(x)$ is defined as
\begin{equation}{h_j}(x) = \frac{{{{( - 1)}^{\kern0.06em j + 1}}(1 - {x^2}){{\tau ^{\prime}}_M}(x)}}{{{M^2}{c_j}(x - {x_j})}},\end{equation}with
\begin{equation}{c_j} = \left\{ {\begin{array}{*{20}{@{}ll}} 1&{1 \le j \le M - 1},\\ 2&{j = 0,M,} \end{array}} \right.\end{equation}
and 
${\tau ^{\prime}_M}(x)$ denotes derivative of the 
$M\textrm{th}$-order Chebyshev polynomial. The interpolant 
${h_j}(x)$ also adheres to the Kronecker delta property:
\begin{equation}{h_j}({x_k}) = {\delta _{jk}},\end{equation}
where 
${\delta _{jk}}$ is known as the Kronecker delta symbol. The explicit determination of derivatives of varying orders can be achieved by differentiating equation (3.3). The expressions for the first- and second-order derivatives at the collocation points 
${x_j}$ are provided below:
\begin{gather}{\left. {\frac{{\textrm{d}f}}{{\textrm{d}\kern0.06em x}}} \right|_j} = \sum\limits_{k = 0}^M {D_{jk}^{(1)}{f_k}} ,\quad j = 0,1, \ldots ,M,\end{gather}
\begin{gather}{\left. {\frac{{{\textrm{d}^2}f}}{{\textrm{d}\kern0.06em {x^2}}}} \right|_j} = \sum\limits_{k = 0}^M {D_{jk}^{(2)}{f_k}} ,\quad j = 0,1, \ldots ,M.\end{gather}
In this context, 
${f_k}$ is essentially just 
$f({x_k})$. The Chebyshev differentiation matrices 
$D_{jk}^{(1)}$ and 
$D_{jk}^{(2)}$ have an order of 
$(M + 1) \times (M + 1)$ and are defined as
\begin{equation}D_{jk}^{(1)} = \left\{
{\begin{array}{*{20}{@{}ll}} {\dfrac{{{c_j}{{( - 1)}^{k +
j}}}}{{{c_k}({x_j} - {x_k})}}}&{j \ne k,}\\ { -
\dfrac{{{x_j}}}{{2(1 - x_j^2)}}}&{1 \le j = k \le M - 1,}\\
{\dfrac{{2{M^2} + 1}}{6}}&{j = k = 0,}\\ { - \dfrac{{2{M^2}
+ 1}}{6}}&{j = k = M,} \end{array}}
\right.\end{equation}
and 
$D_{jk}^{(2)} = D_{jm}^{(1)} \cdot D_{mk}^{(1)}$. The comprehensive methodology is extensively elucidated in the books by Canuto et al. (Reference Canuto, Hussaini, Quarteroni and Zang1988) and Peyret (Reference Peyret2002). The application of this approach to (2.24) and (2.25) gives
\begin{gather}\sum\limits_{k = 0}^M {D_{jk}^{(2)}{\varPsi _k}} - {a^2}{\varPsi _j} - \sum\limits_{k = 0}^M {\sum\limits_{i = 1}^3 {{R_i}D_{jk}^{(1)}{\varPhi _{ik}}} } = 0,\end{gather}
\begin{gather} {-}\frac{{\textrm{i}a}}{2}{\varPsi _j} + \textrm{i}a\left[ {\frac{{({R_1} + {R_2} + {R_3})}}{2}} \right]{x_j}{\varPhi _{ij}} - \frac{1}{{L{e_i}}}\left[ {\sum\limits_{k = 0}^M {D_{jk}^{(2)}{\Phi_{ik}}} - {a^2}{\Phi_{ij}}} \right]\notag\\ = \textrm{i}ac{\varPhi _{ij}}\quad (i = 1,2,3).\end{gather}
Equations (3.10) and (3.11) are now discretized for inner collocation points (
$\,j = 1$ to 
$M$), and the boundary conditions are discretized for 
$j = 0\textrm{ }(x = 1)$ and 
$j = M\textrm{ }(x =- 1)$. This procedure culminates in solving a generalized eigenvalue problem:
\begin{equation}\boldsymbol{\mathsf{AX}} = c\boldsymbol{\mathsf{BX}},\end{equation}
where c and 
$\boldsymbol{\mathsf{X}}$ are the eigenvalue and the eigenvector of field entities, respectively, while 
$\boldsymbol{\mathsf{A}}$ and 
$\boldsymbol{\mathsf{B}}$ are square and complex matrices of order 
$(4M + 4) \times (4M + 4)$. A complex analogue of the QZ algorithm developed by Moler & Stewart (Reference Moler and Stewart1973) is used to solve the aforementioned eigenvalue problem. The coordinates along the neutral stability curves are determined through an iterative secant method, where either the wavenumber a is held constant while varying 
${R_1}$, or vice versa, with 
${R_1}$ fixed while modifying a. This iteration continues until the condition 
${c_i} = 0$ is fulfilled within a predefined error threshold. The eigenvalue c is the eigenvalue with the largest imaginary part for the particular type of disturbance in question. The errors associated with a and 
${R_1}$ at the neutral point are assessed to be below 0.0001 %. The critical values of a and 
${R_1}$ are established through polynomial interpolation based on data points near the nadir of the neutral curve. The critical Darcy–Rayleigh number 
${R_{1c}} = \min (R_{1c}^S,R_{1c}^T)$, where 
$R_{1c}^S$ and 
$R_{1c}^T$ are the critical Darcy–Rayleigh numbers for the stationary and the travelling-wave disturbances, respectively.
4. Results and discussion
The results are presented exclusively for a KCl–NaCl–sucrose aqueous system for which 
${{\kappa _1} = 1.6 \times {10^{ - 5}}\ \textrm{c}{\textrm{m}^2}\ {\textrm{s}^{ - 1}},\ {\kappa _2} = 1.3 \times {10^{ - 5}}\ \textrm{c}{\textrm{m}^2}\ {\textrm{s}^{ - 1}},\ {\kappa _3} = 0.45 \times {10^{ - 5}}\ \textrm{c}{\textrm{m}^2}\ {\textrm{s}^{ - 1}}}$ (i.e. 
$L{e_1} = 1$, 
$L{e_2} = 1.2307$, 
$L{e_3} = 3.5555$) and also for an aqueous solution of heat–KCl–sucrose system for which 
${\kappa _1} = 1.4 \times {10^{ - 3}}\ \textrm{c}{\textrm{m}^2}\ {\textrm{s}^{ - 1}},{\kappa _2} = 1.6 \times {10^{ - 5}}\ \textrm{c}{\textrm{m}^2}\ {\textrm{s}^{ - 1}}, {\kappa _3} = 0.45 \times {10^{ - 5}}\ \textrm{c}{\textrm{m}^2}\ {\textrm{s}^{ - 1}}$ (i.e. 
$L{e_1} = 1$, 
$L{e_2} = 87.5$, 
$L{e_3} = 311.1111$) (Griffiths Reference Griffiths1979). In the former case, it is seen that the Lewis numbers are not far from unity, as is typical of molecular species in molten metals and magmas as well as of laboratory models using salt–sugar solutions. The Darcy–Rayleigh numbers 
${R_2}$ and 
${R_3}$ are allowed to vary, while 
${R_1}$ is taken as a free parameter. The presence of additional diffusing components instils instability on the basic flow and a systematic study of the topology of neutral stability curves has been undertaken by performing linear stability analysis. The neutrally stable configuration is identified by 
${c_i} = 0$ which delimits the boundary between the regions of parametric stability and instability in the 
$(a,{R_1})$ plane.
4.1. Validation of the code
To validate the numerical code, we calculate the critical values of 
${R_1}$, a and c by systematically adjusting the number of collocation points within the sample sets of governing parameters. The findings presented in table 1 indicate that the convergence of the numerical scheme improves as the number of collocation points increases. It is apparent that employing 
$25$ collocation points yields accurate results up to six decimal places for the KCl–NaCl–sucrose system, regardless of the mode of instability. Conversely, for the heat–KCl–sucrose system, achieving satisfactory convergence in the results depends on the mode of instability and necessitates a greater number of collocation points. In particular, when the instability manifests as a stationary mode, a minimum of 
$M = 30$ points is required, whereas travelling-wave instability demands 
$M = 60$ points, as outlined in table 1. Comparing our results under the limiting condition of 
${R_3} = 0$ (i.e. double-diffusive case) with those obtained by Shankar et al. (Reference Shankar, Naveen and Shivakumara2022) in table 2 reveals a high degree of agreement. This rigorous validation process underscores the reliability of the current computational code.
Convergence process of the Chebyshev collocation method for selected values of 
${R_2}$ and 
${R_3}$ for KCl–NaCl–sucrose and heat–KCl–sucrose aqueous systems.
Comparison of results with those of Shankar et al. (Reference Shankar, Naveen and Shivakumara2022) for double-diffusive system 
$({R_3} = 0)$.
4.2. Double-diffusive case
During the course of the investigation, some interesting results were found for double-diffusive fluid systems and they are presented first to get a clear idea about the present problem. In fact, these findings were not disclosed in any of the previous investigations. Figure 2 shows the evolution of neutral stability curves for different values of 
${R_2}$ when 
${R_3} = 0$ for two values of 
$L{e_2} = 1.2307$ (dashed line) and 
$87.5$ (solid line) corresponding to KCl–NaCl and heat–KCl systems, respectively. From the figure, it is observed that the neutral stability curves are closed, comprising either both stationary and travelling-wave modes or only stationary modes, within which the flow is unstable 
$({c_i} > 0)$, while outside of which it defines the parametric conditions of linear stability 
$({c_i} < 0)$. Furthermore, the instability is possible only in a finite range of negative values of 
${R_2}$ depending on the value of 
$L{e_2}$. The salient implication of the closed neutral curve is the requirement of two values of 
${R_1}$ to completely specify the instability of fluid flow. For 
$L{e_2} = 1.2307$, it is evident that the closed neutral stability curve is formed with two segments of consecutively occurring stationary and travelling-wave modes when 
${R_2} =- 1500$ (figure 2a). However, for 
${R_2} =- 700$, the closed neutral curve consists of only single portions of stationary and travelling-wave modes, wherein the latter mode is reduced notably (figure 2b). With further increasing 
${R_2}$, the closed neutral curve progressively transitions to encompass only the stationary mode (figure 2c–f). Furthermore figure 2(a–f) shows the gradual shrinking of the closed neutral curve with increasing 
${R_2}$, which eventually collapses into a single point and completely disappears at 
${R_2} =- 216$. Whereas, for 
$L{e_2} = 87.5$, the closed neutral stability curves include two distinct portions of stationary and travelling-wave modes, a pattern extending across a wider range of 
${R_2}$ (figure 2a–h). It is observed that the unstable region gradually decreases with increasing 
${R_2}$, ultimately vanishing at 
${R_2} =- 0.54$ (figure 2a–j). Interestingly, in the current scenario, the range of 
${R_2}$ within which the instability persists is increased. Besides, the neutral curves of the KCl–NaCl system are encompassed by that of the heat–KCl system with instability occurring at lower wavenumbers in the latter compared with the former. One common feature observed in all these figures is that the onset of instability arises only through the stationary mode in the case of double-diffusive systems.
Evolution of neutral stability curves for different negative values of 
${R_2}$ for double-diffusive system 
$({R_3} = 0)$ when 
$L{e_2} = 87.5$ (solid line) and 
$L{e_2} = 1.2307$ (dashed line). The black and red lines, respectively, denote the stationary and travelling-wave modes and this convention also applies to figures 3–7, 10, 11 and 14–16.
4.3. Triple-diffusive case
As mentioned earlier, two different sets of transport property ratios are considered in discussing the instability of base flow. The results obtained for these two systems are separately discussed.
4.3.1. KCl–NaCl–sucrose system (Le 2, Le 3 ≳ 1)
(a) Topology of the neutral stability curves
The impact of a third diffusing component on the progression of neutral stability curves is illustrated in figures 3 and 4, where different values of 
${R_2}$, spanning both positive and negative domains, are considered when 
${R_3} =- 100,\ L{e_2} = 1.2307$ and 
$L{e_3} = 3.5555$. In contrast to the double-diffusive case, the base flow is found to be unstable for both positive and negative values of 
${R_2}$ and the shape of neutral stability curves is also somewhat unique. For 
${R_2} = 50$, it is seen that the closed neutral curve is formed by two alternative segments of stationary and travelling-wave modes which exhibit local minimum and the least among the two corresponds to the stationary mode (figure 3a). With increasing 
${R_2}$, the closed neutral curve progressively simplifies, eventually consolidating into a single segment accommodating both stationary mode and travelling-wave mode, with the latter mode diminishing (figure 3b) and even ceasing to exist (figure 3c). A noteworthy transformation unfolds at 
${R_2} = 726$, where both modes manifest as separate closed loops, with the travelling-wave mode positioned beneath the stationary one. This alteration denotes a shift in the mode of instability from the stationary to the travelling-wave mode – an unprecedented finding absent in the context of double-diffusive fluid systems. Such a transition mandates the examination of four 
${R_1}$ values to comprehensively analyse the fluid flow stability. The disconnected neutral curve pertaining to the travelling-wave mode expands and merges with the stationary neutral curve, thereby evolving into a unified loop with increasing 
${R_2}$ (figure 3e–g). Moreover, a closed loop of the stationary neutral curve emerges below the current loop signifying the reversion of instability to the stationary mode when 
${R_2}$ reaches the value 
$837$ (figure 3h). This closed stationary neutral curve continues to grow (figure 3i) until it forms a single loop (figure 3j,k). At 
${R_2} = 1353$, it is seen that yet another closed loop of the travelling-wave mode originates on the right-hand side of the pre-existing loop (figure 3l) but it comes together and looms as a single loop with further increase in 
${R_2}$ (figure 3m–o). Throughout this process (figure 3l–o), the stationary mode consistently remains the dominant mode of instability even though the instability region varies. In general, the global minimum decreases with increasing 
${R_2}$. The behaviour of results for negative values of 
${R_2}$ differs significantly from those obtained for positive values and the same is evident from figure 4. In this figure, the closed neutral stability curves exhibit a consistent pattern, featuring pairs of stationary and travelling-wave modes in succession. In contrast to the previous case, the instability persists exclusively through the stationary mode across all considered values of 
${R_2}$. Moreover, as 
${R_2}$ decreases, there is an increase in the instability region and also in the minimum value of 
${R_1}$.
Evolution of neutral stability curves for different positive values of 
${R_2}$ when 
$L{e_2} = 1.2307$, 
$L{e_3} = 3.5555$ and 
${R_3} =- 100$.
Evolution of neutral stability curves for different negative values of 
${R_2}$ when 
$L{e_2} = 1.2307$, 
$L{e_3} = 3.5555$ and 
${R_3} =- 100$.
The neutral curves obtained for 
${R_3} = 100$, using the same values of 
${R_2}$, 
$L{e_2}$ and 
$L{e_3}$ as those considered in figures 3 and 4, are showcased in figures 5 and 6 for some isolation cases. The observations reveal a complete reversal of the neutral curves compared with the findings obtained for 
${R_3} =- 100$. A prominent deduction drawn from these visualizations is the unswerving nature of the onset of mode of instability in this particular context, which consistently assumes the form of the stationary mode. This is because, regardless of the values of 
${R_2}$, when 
${R_3} =- 100$, the maximum value of 
${R_1}$ reached on the closed neutral curve corresponds to the stationary mode.
Evolution of neutral stability curves for different negative values of 
${R_2}$ when 
$L{e_2} = 1.2307$, 
$L{e_3} = 3.5555$ and 
${R_3} = 100$.
Evolution of neutral stability curves for different positive values of 
${R_2}$ when 
$L{e_2} = 1.2307$, 
$L{e_3} = 3.5555$ and 
${R_3} = 100$.
Figure 7 presents a comprehensive depiction of the neutral stability curves for different positive values of 
${R_2}$ when 
${R_3} =- 200$. The neutral curve loop unravels intricate dynamics involving both stationary and travelling-wave modes for 
${R_2} = 1050$, featuring two minima wherein the stationary mode predominates as the mode of instability (figure 7a). Evidencing further complexity, at 
${R_2} = 1110$, the neutral curve portrays three local minima, with the initial minimum corresponding to the stationary mode and the subsequent two to the travelling-wave mode. Notwithstanding this complex scenario, the mode of instability continues as the stationary mode (figure 7b). Intriguingly, as 
${R_2}$ increases to 
$1150$, a transition emerges, wherein the mode of instability transitions from the stationary mode to the travelling-wave mode (figure 7c). A noteworthy change in the nature of the neutral stability curves manifests at 
${R_2} = 1303$, as they give up their multi-minimum structure, rendering a single minimum that aligns with the travelling-wave mode (figure 7d). Of particular interest is the transition observed at 
${R_2} = 1304$, marked by the emergence of an extra closed loop depicting the stationary mode (figure 7e), necessitating the consideration of four critical values of the Darcy–Rayleigh number to assess the stability of the system. Significantly, the stationary mode reasserts itself as the favoured mode of instability and the closed loop of the stationary mode grows with further increase in 
${R_2}$ (figure 7f). Importantly, within the intricate interplay of these loops, another closed loop emerges, characterizing the travelling-wave mode (figure 7g). This emergent loop undergoes gradual evolution and subsequently establishes a connection with the lower stationary loop (figure 7h,i), leading to their integration into a coherent single loop (figure 7j) as 
${R_2}$ progresses. Moreover, yet another closed loop representing the travelling-wave mode appears independently without altering the mode of instability (figure 7k,l) and ultimately amalgamates with the pre-existing closed loop (figure 7m–o). This cumulative loop entails multiple segments, encompassing both stationary and travelling-wave modes. The gradual transformation of these loop structures reveals the complex behaviour of the system in response to varying 
${R_2}$ values, highlighting the interplay between different modes of instability.
Evolution of neutral stability curves for different positive values of 
${R_2}$ when 
$L{e_2} = 1.2307$, 
$L{e_3} = 3.5555$ and 
${R_3} =- 200$.
(b) Stability boundaries
The variation of critical Darcy–Rayleigh number 
${R_{1c}}$ as a function of 
${R_2}$ for different values of 
${R_3}$ when 
$L{e_2} = 1.2307$ and 
$L{e_3} = 3.5555$ is displayed in figure 8. It is clear that the effect of increasing both 
${R_2}$ and 
${R_3}$ is to decrease 
${R_{1c}}$ indicating their effect is to destabilize the fluid flow. For the double-diffusive case 
$({R_3} = 0)$, instability is solely possible through the stationary mode up to 
${R_2} =- 216.08$ and thereafter the flow remains stable. In contrast to this, it is observed that the instability exists for all values of 
${R_2}$ if the third diffusing component is present. Moreover, the transition of instability changes from the stationary mode (denoted by solid lines) to the travelling-wave mode (denoted by dotted lines) and again reverts back to the stationary mode with increasing 
${R_2}$, a phenomenon intricately linked to the value of 
${R_3}$. In particular, the travelling-wave mode instability is present only for negative values of 
${R_3}$ and it is confined to a finite domain of positive values of 
${R_2}$ as illustrated in figure 8(b). Also note that this range of 
${R_2}$ initially increases with decreasing 
${R_3}$ then starts decreasing and eventually converges to zero as 
${R_3}$ experiences further reduction. The inherent instability characteristics of fluid flow, manifested across distinct ranges of the parameter 
${R_2}$, corresponding to different values of 
${R_3}$, have been discretely tabulated in table 3.
Critical values of 
${R_1}$ versus 
${R_2}$ for both positive (a) and negative (b) values of 
${R_3}$ when 
$L{e_2} = 1.2307$ and 
$L{e_3} = 3.5555$. The solid line represents a stationary mode, while the dotted line corresponds to a travelling-wave mode and this convention is also applicable to figures 9, 12 and 13.
Effect of 
${R_3}$ on the instability transitions from one mode to another in different ranges of 
${R_2}$ for the KCl–NaCl–sucrose aqueous system.
The variation of corresponding critical wavenumber 
${a_c}$ with 
${R_2}$ is presented in figure 9. For positive and also for higher negative values of 
${R_3}$, no discontinuity in the curves of 
${a_c}$ exists as the instability occurs only through the stationary mode in the considered range of values of 
${R_2}$. In particular, for 
${R_3} = 100$ and 
$200$, 
${a_c}$ passes through a maximum value with increasing 
${R_2}$, while for other positive values of 
${R_3}\textrm{ }( = 300,400,1000)$, it decreases for all values of 
${R_2}$ but the trend is found to be not so significant. An entirely distinct scenario becomes visible for negative values of 
${R_3}\textrm{ }( =- 100, - 200, - 300, - 400)$, a phenomenon attributed to the presence of transition modes, which in turn leads to the emergence of discontinuities in 
${a_c}$. For 
${R_3} =- 1000$, however, no transition occurs and there is only a slight increase in 
${a_c}$ throughout the range of considered values of 
${R_2}$. The critical wavenumber corresponding to the stationary mode increases prior to the initiation of the travelling-wave mode, while in post-transition, this 
${a_c}$ value diminishes. It is pertinent to note that within the travelling-wave mode regime, 
${a_c}$ registers a decreasing trend. In the framework of the double-diffusive context 
$({R_3} = 0)$, 
${a_c}$ increases with increasing 
${R_2}$, reaches the maximum value and then a downturn in 
${a_c}$ follows as 
${R_2}$ further increases and finally terminates at 
${R_2} =- 216.08$. The critical wave speed 
${c_c}$ increases as 
${R_2}$ increases; concurrently with a decrease in 
${R_3}$ values (figures are not shown). Notably, travelling-wave disturbances propagate both upward and downward vertically at the same velocity, owing to the existence of both positive and negative 
${c_c}$ values.
Critical values of a versus 
${R_2}$ for both positive (a) and negative (b) values of 
${R_3}$ when 
$L{e_2} = 1.2307$ and 
$L{e_3} = 3.5555$.
4.3.2. Heat–KCl–sucrose system 
$(L{e_2},\ L{e_3} \gg 1)$
(a) Topology of the neutral stability curves
The neutral stability curves are shown in figure 10 for different positive values of 
${R_2}$ when 
${R_3} =- 1000,\ L{e_2} = 87.5$ and 
$L{e_3} = 311.1111$. As observed in the previous case, the neutral curves form closed and connected patterns, consisting of both stationary and travelling-wave modes. For 
${R_2} = 1000$, multiple minima linked to both modes are evident, with primary minimum corresponding to the stationary mode, which is identified as the favoured mode of instability (figure 10a). Interestingly, two segments of the neutral stability curve of stationary mode possessing their own local minimum appear for 
${R_2} = 1050$ and again the instability occurs at the initial wavenumber region (figure 10b). Analogously, when 
${R_2} = 1100$, the neutral stability curves follow the same pattern (figure 10c), but the primary instability is now dominated by the local minimum occurring at a higher wavenumber of the stationary mode neutral curve. This particular mode maintains its dominance as the preferred mode of instability even with increased values of 
${R_2}$ (figure 10d,e).
Evolution of neutral stability curves for different positive values of 
${R_2}$ with 
$L{e_2} = 87.5,\textrm{ }L{e_3} = 311.1111$ and 
${R_3} =- 1000$.
Figure 11 portrays the neutral stability curves for different negative values of 
${R_2}$ when 
${R_3} = 200,\ L{e_2} = 87.5$ and 
$L{e_3} = 311.1111$. We note that a single loop of neutral stability curve featuring multi-minima associated with both modes can be seen with the stationary mode as the preferred mode of instability for 
${R_2} =- 344$ (figure 11a). However, an entirely detached closed travelling-wave neutral curve emanates for 
${R_2} =- 343$ which lacks any physical significance for the onset of instability of the system as the minimum 
${R_1}$ still corresponds to the stationary mode (figure 11b). In the subsequent cases, the detached travelling-wave loop gradually expands (figure 11c) and eventually merges (figure 11d–f) with the pre-existing neutral curve with increasing 
${R_2}$. Remarkably, during this process, a significant shift takes place in the instability mode, transitioning from the stationary mode to the travelling-wave mode (figure 11e). However, as 
${R_2}$ continues to increase, the system reverts back favouring the stationary mode of instability once again (figure 11f).
Evolution of neutral stability curves for different negative values of 
${R_2}$ with 
$L{e_2} = 87.5,\textrm{ }L{e_3} = 311.1111$ and 
${R_3} = 200$.
(b) Stability boundaries
Figures 12 and 13 demonstrate, respectively, the variation of 
${R_{1c}}$ and 
${a_c}$ as a function of 
${R_2}$ for different values of 
${R_3}$ when 
$L{e_2} = 87.5$ and 
$L{e_3} = 311.1111$. In the figures, the curve of 
${R_3} = 0$ corresponds to the double-diffusive system and it is observed that 
${R_{1c}}$ tends to infinity for values of 
${R_2} > - 0.541$ indicating that the flow is always stable, because the perturbations display a negative growth rate. While for values of 
${R_3} \ne 0$, the instability is found to pervade the entire 
${R_2}$ domain, as depicted in both panels of figure 12. More specifically, the instability predominantly manifests through the travelling-wave mode (illustrated by dotted lines) within a specific interval of negative 
${R_2}$ for positive 
${R_3}$ values, while in other regimes, instability exclusively arises through the stationary mode (illustrated by solid lines). Additionally, it is noteworthy that the span of 
${R_2}$ values, where instability is induced through travelling waves, diminishes as 
${R_3}$ decreases. Moreover, it is worth noting that increasing both 
${R_2}$ and 
${R_3}$ serves to speed up the onset of thermoconvective instability in both stationary and travelling-wave modes. The salient observations discussed above are quantified in table 4, providing a comprehensive summary. The critical wavenumber 
${a_c}$ demonstrates distinct trends with increasing values of either 
${R_2}$ or 
${R_3}$ (figure 13). A noteworthy observation for positive values of 
${R_3}$ is the abrupt jump in 
${a_c}$ values within the negative domain of 
${R_2}$, signifying a pivotal transition in the instability mode from stationary to travelling wave and this behaviour is confirmed quite clearly by the neutral curves displayed in figure 11. Also, during the travelling-wave mode, the variation in 
${a_c}$ becomes relatively less pronounced as 
${R_2}$ increases. While for negative values of 
${R_3}$, a sudden variation in 
${a_c}$ curves ensues in the positive domain of 
${R_2}$ due to the shift in the mode of stationary instability from lower- to higher-wavenumber region (cf. figure 10). The critical wave speed 
${c_c}$ increases proportionally with increased values of 
${R_2}$ and 
${R_3}$ (figures are not shown).
Critical values of 
${R_1}$ versus 
${R_2}$ for both positive (a) and negative (b) values of 
${R_3}$ when 
$L{e_2} = 87.5$ and 
$L{e_3} = 311.1111$.
Critical values of a versus 
${R_2}$ for both positive (a) and negative (b) values of 
${R_3}$ when 
$L{e_2} = 87.5$ and 
$L{e_3} = 311.1111$.
Effect of 
${R_3}$ on the instability transitions from one mode to another in different ranges of 
${R_2}$ for the heat–KCl–sucrose aqueous system.
4.3.3. Comparative study between KCl–NaCl–sucrose and heat–KCl–sucrose systems
A comprehensive comparative analysis between the KCl–NaCl–sucrose and heat–KCl– sucrose systems has been undertaken to thoroughly understand both the shared characteristics and unique distinctions for the stability of fluid flow. This is facilitated by presenting neutral stability curves in figures 14–16, for identical values of 
${R_2}$ and 
${R_3}$. Within these graphical representations, the dashed and solid lines symbolize the KCl–NaCl–sucrose and heat–KCl–sucrose systems, respectively. The neutral curves of the KCl–NaCl–sucrose system either fall comfortably within or intersect with those of the heat–KCl–sucrose system. However, the instability region delineated by the neutral curves of the KCl–NaCl–sucrose system is significantly smaller when compared with the latter system. Additionally, it is observed that the onset of instability happens first in the heat–KCl–sucrose system, no matter what values the governing parameters take. A keen scrutiny of the aforementioned figures also indicates that the transition mode of instability for these two aqueous systems takes place in different regimes of 
${R_2}$ and 
${R_3}$. Furthermore, an observation of significance emerges, indicating that the KCl–NaCl–sucrose system necessitates the identification of four critical values of the Darcy–Rayleigh number to establish the linear stability criteria under specific parametric conditions. In stark contrast, the heat–KCl–sucrose system mandates only two critical values for this purpose, irrespective of the values assumed by 
${R_2}$ and 
${R_3}$.
Comparison of neutral stability curves between KCl–NaCl–sucrose (dashed line) and heat–KCl–sucrose (solid line) aqueous systems for different positive values of 
${R_2}$ with 
${R_3} =- 100$.
Comparison of neutral stability curves between KCl–NaCl–sucrose (dashed line) and heat–KCl–sucrose (solid line) aqueous systems for different positive values of 
${R_2}$ with 
${R_3} =- 1000$.
Comparison of neutral stability curves between KCl–NaCl–sucrose (dashed line) and heat–KCl–sucrose (solid line) aqueous systems for different negative values of 
${R_2}$ with 
${R_3} = 200$.
5. Conclusions
The linear stability of natural convection in a vertical triple-diffusive fluid-saturated porous layer is investigated. The dynamics of convection is governed by Darcy's law, while the buoyancy effects are characterized by employing the Oberbeck–Boussinesq approximation. A numerical solution of the stability eigenvalue problem has been obtained based on the Chebyshev collocation method for two different transport property ratios relating to KCl–NaCl–sucrose and heat–KCl–sucrose aqueous systems. The evolution of neutral stability curves is examined in detail and it is observed that these two systems exhibit qualitatively different behaviour for the stability of base flow under different regimes of solute Darcy–Rayleigh numbers 
${R_2}$ and 
${R_3}$. Besides, the presence of a third diffusing component not only renders the system unstable without any restriction, in contrast to double-diffusive systems, but it also introduces a transition mode that is found to be absent within the ambit of double-diffusive systems.
Some of the novel peculiar features gathered from the foregoing study under different regimes of 
${R_2}$ and 
${R_3}$ may be summarized as follows:
(i) When
${R_2}$ and ${R_3}$ are concurrently positive or negative, the neutral stability curves exhibit a single loop encompassing stationary and travelling-wave modes and the instability occurs only through the stationary mode for both the aqueous systems (i.e. no transition mode occurs). Also, two critical values of the Darcy–Rayleigh number are needed to dwell upon the stability criteria. This phenomenon is akin to that observed in double-diffusive fluid systems.(ii) For values of
${R_2} < 0$ and ${R_3} > 0$, the possibility of two entirely disjointed neutral curves emerging is found for the KCl–NaCl–sucrose system indicating the necessity of four critical values of the Darcy–Rayleigh number to precisely delineate the stability criteria and the instability exclusively materializes in a stationary mode. On the contrary, the heat–KCl–sucrose system exhibits the sufficiency of two critical values of the Darcy–Rayleigh number for the precise specification of stability criteria despite the neutral curves demonstrating multiple minima and going through a transition mode. The range of ${R_2}$ at which the travelling-mode occurs increases with increasing ${R_3}$.(iii) When
${R_2}$ takes on positive values and ${R_3}$ adopts negative values, the neutral curves are totally disconnected under a certain parametric condition as in (ii) for the KCl–NaCl–sucrose system. Yet, this shift in configuration precipitates a transition in the mode of instability. Specifically, a regime of travelling-wave mode emerges within a defined domain of positive ${R_2}$ which progressively diminishes as ${R_3}$ decreases further. While for the heat–KCl–sucrose system, only a single loop structure emerges without any transition mode.
Thus the presence of additional diffusing components and the magnitude of transport property ratios assume a pivotal role in delineating the stability characteristics of the system. The instability region outlined by the neutral curves of the KCl–NaCl–sucrose system is notably smaller in contrast to the heat–KCl–sucrose system. Furthermore, the heat–KCl–sucrose system demonstrates greater destabilization compared with the KCl–NaCl–sucrose system. A common feature across both the systems is that increasing the values of either 
${R_2}$ or 
${R_3}$ advances the onset of instability, regardless of the mode of instability.
The findings discussed in this study disclose a spectrum of potential avenues for future advancements. One possibility of exploration entails the conceptualization of a sandwiched porous slab, characterized by uniform external layers but varying thermophysical properties within the core layer. This aspect holds significant importance as it offers a deeper insight into convection dynamics in complex configurations. Additionally, a separate research initiative will aim to extend the current work by incorporating the general boundary conditions on velocity, temperature and/or solute concentrations. Another notable concern revolves around the nonlinear stability which may reveal the potential emergence of subcritical instability, providing a valuable comparison to establish the stability threshold.
Acknowledgements
We express our gratitude to the Associate Editor and an anonymous referee for their valuable suggestions, which led to improvements in the manuscript.
Declaration of interests
The authors report no conflict of interest.
Data availability
The data that support the findings of this study are available within the article.
