Analytical and Numerical Study of Soret and Dufour Effects on Thermosolutal Convection in a Horizontal Brinkman Porous Layer with a Stress-Free Upper Boundary

In this paper, thermo-diffusion (Soret effect) and diffusion-thermo (Dufour effect) effects on double-diffusive natural convection induced in a horizontal Brinkman porous layer with a stress-free upper boundary are investigated.)e cavity is filled with a binary fluid and subjected to uniform fluxes of heat and mass on its long sides. An analytical solution based on the parallel flow approximation is developed for the problem considered in order to allow prompt determination of the thresholds of stationary and finite amplitude solutions and also heat and mass transfer characteristics. )e analytical solution is validated numerically by using a finite difference method. )e combined effects of the Soret and Dufour parameters, the thermal Rayleigh number, the buoyancy ratio, and the Darcy number on the flow intensity and heat and mass transfer are illustrated graphically, and some particular behaviors observed are discussed.)e analytical solution proves the existence of different regions in the buoyancy ratioDufour parameter plane, corresponding to different parallel flow behaviors. )e number, the location, and the extent of these regions, which are impossible to predict numerically, depend strongly on Soret and Dufour parameters. )e effect of thermodiffusion and diffusion-thermo on flow intensity and heat and mass transfer is found to be important.


Introduction
Natural convection involving binary mixtures in the porous and fluid media is still a relevant field of investigation since it occurs in wide ranges of applications in many engineering problems and natural fields such as geophysics, oil reservoirs, storage of nuclear wastes, operation of solar ponds, chemical reactors, migration of mixtures in fibrous insulation, and metal manufacturing processes. is phenomenon occurs due to temperature differences, concentration differences, or by combination of these two differences [1]. Convection in which the buoyancy forces are due to both temperature and chemical concentration gradients is referred to as thermosolutal or double-diffusive convection. e fluid flow behavior driven by temperature and concentration gradients has attracted the interest of researchers worldwide through the decades [2]. Different aspects of the problems involving double-diffusive natural convection have been addressed previously [3][4][5][6][7][8][9][10] in the absence of Soret and Dufour effects and neglected also in many other studies since they are of smaller order of magnitude compared to the effects described by Fourier's and Fick's laws. Nevertheless, the literature review shows the existence of experiments mainly focused on the determination of the Soret coefficients for mixtures [11][12][13] or on theoretical aspects to determine the optimal combinations of the governing parameters leading to the separation of species under the effect of thermo-diffusion [14][15][16]. Other aspects of double-diffusive convection problems in the presence of the Soret effect have been conducted in shallow horizontal cavities by Bourich et al. [17][18][19]. e results reported in these investigations are concerned with convection thresholds [17], reversal of the horizontal concentration gradient [18], and a comparison study between the cases where the enclosure is filled with a clear binary mixture or a saturated porous medium [19]. More recently, Hasnaoui et al. [20] studied numerically double-diffusive natural convection in an inclined enclosure with heat generation and Soret effect using a hybrid lattice Boltzmann-finite difference method. ey show that the negative Soret parameter combined with high internal heat generation and a relatively high inclination is important when the objective is to maintain the fluid at a high concentration of species.
Despite the fact that Dufour and Soret effects are of a smaller order of magnitude, they however may play a significant role in double-diffusive flow processes. e list of published works taking into account both these two effects is limited. Ojjela et al. [21] studied numerically the influence of an external magnetic field and radiation in free convective Jeffrey fluid flow between parallel porous plates in the presence of Soret and Dufour effects. e governing equations were solved numerically using the shooting method with the 4 th Runge-Kutta scheme. e numerical study conducted by Zaho et al. [22] focused on heat and mass transfer generated by natural convection in a porous medium saturated by a viscoelastic fluid and subjected to a magnetic field. Kefayati [23] studied the impact of Soret and Dufour effects on entropy generation due to double diffusion in a square cavity filled with non-Newtonian powerlaw fluid by using a finite-difference lattice Boltzmann method. He found that the Dufour parameter results in the enhancement of the total irreversibility. A similar study was performed by the same author [24,25] in the case of inclined cavities. Wang et al. [26] investigated numerically the behavior change of solutions accompanying the increase of the buoyancy ratio. Ren and Chan [27] studied transient double-diffusive convection under the effect of numerous parameters of control including thermodiffusion and diffusion-thermo. Balla and Naikoti [28] analyzed numerically double-diffusive free convection heat and solute transfer in an inclined square porous cavity saturated with a fluid in the presence of Soret and Dufour effects. Both these effects were considered by Al-Mudhaf et al. [29] who studied numerically unsteady double-diffusive natural convection inside trapezoidal inclined enclosures filled with an isotropic porous medium and submitted from one side to nonconstant distributions of temperature and concentration. Double-diffusive natural convection of non-Newtonian power-law fluids in an open cavity subjected to a horizontal magnetic field in the presence of Soret and Dufour effects was analyzed by Kefayati [30,31]. Among others, he showed that the rise of Soret and Dufour parameters enhances the entropy generations due to heat transfer and fluid friction. e present study is devoted to analytical and numerical investigation of Soret and Dufour effects within a horizontal Brinkman layer with a stress-free upper boundary. e main objective is to bring out the combined effects of Soret and Dufour parameters on thresholds of convection, fluid flow, and heat and mass transfer characteristics. For this shallow enclosure, the parallel flow approximation allows to predict the critical Rayleigh numbers for the onset of supercritical and subcritical convection. In addition, the simplifications allowed by this approximation leads to identify different regions with various parallel flow behaviors. e adopted Brinkman model permits covering the limits of Darcy and the pure fluid media. e study focuses on the discussion of combined effects of Rayleigh number, Soret parameter, Dufour parameter, buoyancy ratio, and Darcy number.

Mathematical Formulation of the Problem
e domain under study, sketched in Figure 1, is a two-dimensional horizontal porous cavity saturated with a binary mixture. e height of the cavity is H ′ , its length is L ′ (A � L ′ /H ′ ≫ 1), and its upper horizontal surface is nondeformable and stress-free, while the remaining boundaries are assumed rigid. e short boundaries (vertical walls) of the cavity are maintained adiabatic and impermeable to mass transfer, while its long horizontal walls are subjected to uniform fluxes of heat, q ′ , and mass, j ′ . e porous matrix is assumed isotropic and homogeneous, and the Brinkman-Hazen-Darcy model is adopted. e binary fluid that saturates the porous matrix is assumed to be Newtonian, and it is modeled as a Boussinesq incompressible fluid whose density ρ varies linearly with the temperature T ′ and the concentration S ′ as follows: where ρ 0 is the binary fluid density at temperature and concentration (T 0 ′ , S 0 ′ ) and β T and β S are the thermal and solutal expansion coefficients, respectively. e dimensional equations governing this problem are written as e reference scales for time, length, velocity, and pressure are σH ′ 2 /α, H ′ , α/H ′ , and αμ/K, respectively. e parameters σ and α are defined by σ � (ρC) P /(ρC) f and α � λ/(ρC) f . e dimensionless temperature and mass fraction are e governing equations that are solved numerically are based on the vorticity-stream function formulation. In their dimensionless forms [26], they are obtained as follows: 2 Mathematical Problems in Engineering e boundary conditions associated to the present problem are e examination of the dimensionless governing equations, equations (6)-(10), and their associated boundary conditions, equations 11a-11c, shows that the control parameters are the aspect ratio of the enclosure, A, Darcy number, Da, the Dufour number, Du, the Lewis number, Le, the thermal Rayleigh number, R T , and the Soret number, Sr. ese parameters are defined as follows: e Nusselt and Sherwood numbers are, respectively, given by the following expressions: where ΔT � T(0, 1/2) − T(0, − 1/2) and ΔS � S(0, 1/2) − S (0, − 1/2) are the temperature and concentration differences, evaluated at x � 0.

Numerical Solution
e numerical solution of the governing equations together with the associated boundary conditions, equations (6) to (11a)-(11c), is obtained using a finite-difference method. e vorticity equation, equation (6), the energy equation, equation (7), and the equation of conservation of species, equation (8), are written in their transient forms since the iterative procedure is performed using the alternate direction implicit (ADI) method for these equations. e stream function field is derived from the resolution of equation (9) using the point successive over-relaxation method (PSOR). e numerical results reported in this paper were obtained with a globally nonuniform grid; the latter is uniform and tight near the confining rigid walls and the free surface to capture the flow details near these boundaries and uniform but coarser elsewhere. e vorticity values on the rigid boundaries are calculated using Wood's relation [32], while zero is attributed to the vorticity at the free surface. Further details on the numerical method and its validation were given in a previous study by Amahmid et al. [5].
Numerous preliminary numerical tests have been performed to determine the minimum aspect ratio from which the assumption of the parallel flow is corroborated. Within the ranges of variations of the parameters considered in this investigation, it was found that the numerical results are independent of the aspect ratio from the threshold A � 12. Mathematical Problems in Engineering ereby, the numerical results reported here were obtained with A � 12 and a grid of 201 × 81. e effect of the grid size on the results is illustrated in Table 1 (analytical results are also included in this table). It can be seen from this table that the results obtained with the grid 201 × 81 differ by less than 0.2% from those corresponding to the finest grid 261 × 121. In addition, the grid 201 × 81 reproduces the analytical results with relative errors lower than 0.3%.

Parallel Flow Analysis.
Based on the remarks (verified numerically) allowing simplifications in the case of shallow cavities (A ≫ 1), the flow is parallel to the long walls of the cavity and the temperature and concentration fields are characterized by a linear and horizontal stratification if we exclude the edges' effects.
us, the stream function, the temperature, and the concentration can be approximated as follows: the parameters C T and C S are the unknown constant temperature and concentration gradients in the horizontal direction, respectively. e investigation by Cormack et al. [33] counts among the early works developing these approximations based on the theory of asymptotic developments for A ≫ 1. e steady form of the governing equations (6)-(8) while using these approximations can be reduced to the following set of ordinary differential equations: where R � − R T (C T + NC S ). e corresponding boundary conditions on the horizontal walls become e analytical integration of equations (15)- (17), together with the associated boundary conditions (18) and (19), leads to the following solutions, respectively, in terms of stream function, horizontal velocity, temperature, and concentration: e parameter ψ 0 is the stream function at the center of the cavity, obtained analytically as e parameters in the previous equation are obtained as is the normalized Rayleigh number, and Z(y), H(y), and S(y) are the functions, respectively, given by the following analytical expressions: e expressions of χ i (i � 1, 2, 3, 4) are given by e parameter ζ � Da − 1/2 and the functions A(Da) and Γ(Da) are dependent on the Darcy number as follows: According to equation (21), the velocity cannot be zero in the range 0 < y < 1/2, which means that only unicellular flows are possible. In addition, the constants C T and C S are obtained by considering mass and thermal balances across any transversal section of the porous layer [34].
Substituting equations (22) and (23) into equation (13), the respective analytical expressions of Nusselt and Sherwood numbers are obtained as By substituting the solutions given by equations (21)-(23) into equations (31) and (32), the final analytical expressions of C T and C S are obtained as follows: with M � 1 − Du Sr and R � 2ψ 0 /Z(0).
Mathematical Problems in Engineering e analysis of equation (24) indicates that up to five different steady-state solutions are possible including the diffusive regime (ψ 0 � 0). e four other convective solutions depend on the signs ± inside and outside the brackets. In fact, the sign ± outside the brackets refers to counterclockwise/clockwise circulation, while it refers to convective stable/unstable solutions within the brackets. It is to specify that the "unstable" solutions refer to the analytical solutions that could not be obtained numerically. In addition, from a mathematical point of view, equation (24) shows the existence of supercritical and subcritical bifurcations [19]. e supercritical bifurcation is characterized by the transition from the quiescent state to the convective regime through zero flow amplitude (ψ 0 � 0), while the subcritical bifurcation is characterized by the onset of motion through finite amplitude ψ 0 � ± Z(0) . ese parallel flow solutions exist only when the two following conditions are satisfied: e satisfaction of these conditions allows the subdivision of the (N, Du) plane into different regions with specific behaviors, depending on the sign of the parameter P defined by equation (26).  (1) is characterized by the rest state, which means that the parallel flow solution is not possible in this region whatever the value of R T . Different conditions must be verified to delineate this region in the (N, Du) plane, depending on the sign of Sr.
For the region (2), only the supercritical flow is possible. e corresponding supercritical Rayleigh number, characterizing the transition from the rest state to the convective regime through the zero flow amplitude, is given by Region (1) Region (5) Region (3) Region (2) Region (4)  Region (5) Region (4) Region (1) Region (2) Region (3) Mathematical Problems in Engineering e conditions that should be verified to delineate this region in the (N, Du) plane are dependent on the sign of Sr.
For the region (3), the characteristics are similar to those in region (2). e stable parallel flow solution exists from the Region (4) Region (5) Region (1) Region (2) Region ( (2) and (3) correspond to the same type of convection (stationary convection), they are however characterized by a difference in terms of asymptotic evolutions at large values of the Rayleigh number. e conditions that should be verified for this region are dependent on the sign of Sr as follows: and , For the region (4), both stable and unstable solutions are existing from the subcritical Rayleigh number, but the unstable branch disappears at some threshold of R T . In this region, the resulting solutions correspond to a subcritical bifurcation for which the onset of motion occurs through a finite amplitude. e critical Rayleigh number, R sub TC , above which the parallel flow exists, is obtained from the condition b 2 − 4Le 2 c � 0 and its expression is given by where G � (Le 4 − Le 2 )+ N(Le − Le 3 ) + 2N Sr Du (Le 3 − Le 2 ) + (2N Sr + 2Du)(Le 2 + Le 3 ). e conditions that should be verified for the region (4) in the (N, Du) plane are dependent on the sign of Sr as follows: 8 Mathematical Problems in Engineering For the region (5), the behaviors are similar to those of region (4) but the existence range of the unstable branches of the former is extended towards infinite R T .
Finally, the region (6) is characterized by the same behavior as that of region (4); the only difference is observed in terms of their asymptotic behaviors.
e Soret number appears explicitly in the conditions delineating this region as follows: Mathematical Problems in Engineering e expression of f(Du) was already defined for region (2).   corresponding to different regions in this plane. e choice of Sr � 0.6 allows to cover the six identified regions. Moreover, it can be observed from these figures that the stable analytical solution (solid lines) is in a very good agreement with the numerical results depicted by full circles and obtained by solving the full governing equations. For the combination (N, Du) � (− 2, − 0.5) corresponding to the region (2), only the stable branch exists. e supercritical convection starts from the rest state (characterized by ψ 0 � ψ(0, 0) � 0, Nu � 1, and Sh � 1) at R Sup TC � 11.12. e evolutions of ψ 0 , Nu, and Sh are characterized by asymptotic behaviors, and their respective asymptotic limits are 2.31, 3.45, and 1.64, reached for sufficiently large values of R T . Note that for this combination of (N, Du), the asymptotic value of Sh indicates that diffusion plays an important role in the mass transfer. At large values of R T , the analytical expressions of the asymptotic expressions of ψ 0 , Nu, and Sh corresponding to the stable solution are obtained by the following expressions:

Effect of
For the combination (N, Du) � (2, 0.7), illustrating the region (3), the behavior is also supercritical and the convection starts from the rest state at R Sup TC � 11.42. e evolution of ψ 0 vs. R T is characterized by a monotonous increase with different rates that depend on the combination (N, Du). For R T ⟶ + ∞, the analytical expression of ψ 0 is reduced to Equation (58) shows that ψ 0 varies as R 1/2 T , while the evolutions of Nu and Sh are characterized by asymptotic behaviors towards the same limit, which is 4.49 (see Figures 5(b) and 5(c)). is value is deduced from the asymptotic expression given by equation (59) that is common for both Nu and Sh at large R T : (59) e region (4) is illustrated by the combination (N, Du) � (0.3, 0.7) in Figures 5(a)-5(c). Unlike the previous combinations of (N, Du), for which only the stable branches exist, both stable and unstable solutions exist for this combination in the range R Sub TC � 45.12 ≤ R T < R Sup TC � 83.50. is means that the parallel flow solution starts at R Sub TC � 45.12, and the nascent flow corresponds to a state clearly different from the purely diffusive regime characterized by ψ 0 � 0 and Nu � Sh � 1. For the unstable branch, ψ 0 /(Nu and Sh) decreases notably/(weakly) by increasing R T and disappears for R T ≥ R Sup TC � 83.50.
For the stable branch, the evolution of ψ 0 versus R T is very similar to that described for the combination (N, Du) � (2, 0.7) but characterized by a delay in terms of R T . At large values of the last parameter, the expression of ψ 0 is the same as that established for the region (3) and given by equation (58). For the region (4) also, the evolutions of Nu and Sh are characterized by asymptotic behaviors at large R T , leading to the common asymptotic value 4.49, obtained with equation (59) for both Nu and Sh.
e results corresponding to the region (5) are illustrated in Figures 5(a)-5(c) by the combination (N, Du) � (0.2, 0.9). For this region, both stable and unstable solutions exist for R T ≥ R Sub TC � 69.16 for the considered combination of (N, Du). In addition, the evolution of the stable branch is similar to that described for region (4). By increasing R T , the unstable branches are characterized by slow decreases towards asymptotic limits for ψ 0 , Nu, and Sh. e region (6) is illustrated by the combination (N, Du) � (− 1.2, − 1.8) and shows a behavior different from those already described. More precisely, both stable and unstable solutions for this region exist in the very short range R Sub TC � 14.793 ≤ R T < R Sup TC � 14.797. In addition, at large values of R T , the quantities ψ 0 , Nu, and Sh tend, respectively, to the asymptotic limits 9.70, 5.827, and 4.155. It should be outlined that for the region (6) the behaviors are similar to those already described for the regions (4) and (2), respectively, at low and large values of R T . Consequently, the analytical expressions of the asymptotic values of ψ 0 , Nu, and Sh are given by equations (55)-(57), respectively. Besides, for the region (6), the Nusselt number becomes negative in the range 27.34 ≤ R T ≤ 453.99. is change in the sign is exclusively attributed to the negative value of the Dufour number and the small temperature differences induced in some ranges of R T .

Combined Effects of Dufour and Soret Numbers.
is section is devoted to assess the simultaneous effects of thermal-diffusion (Soret effect) and diffusion-thermo (Dufour effect) on fluid flow and heat and mass transfer characteristics. Soret and Dufour effects are defined as the mass flux caused by a temperature difference and the energy flux engendered by concentration difference, respectively [35]. e Soret and Dufour parameters can be widely varied by changing the mean solute concentration or mean temperature of the system. e effect of the Dufour number on the flow intensity, the Nusselt number, and the Sherwood number is illustrated e effect of Sr on ψ 0 , Nu, and Sh may be drastically affected by the change of Du. More specifically, we can observe two different behaviors depending on whether the value of Du is above or below critical ranges, which are so narrow that they appear as nodes in Figures 6(a)-6(c). In these critical ranges of Du, the effect of thermo-diffusion becomes insignificant. Another mystery of these critical ranges of Du lies in the fact that they are not the same for ψ 0 , Nu, and Sh (compare [0.6, 0.732] for ψ 0 , [0.016, 0.042] for Nu, and [− 2.29, − 2.25] for Sh). Due to their very narrowness, the critical ranges will be termed as critical nodes in the following discussion. e critical node corresponding to the curves of Sh is seen to be largely shifted towards negative values of Du. e interaction between convection, themodiffusion, and diffusion-thermo is so complex that, without analytical calculations, it was not possible to predict such specific behaviors attributed to the combined effects of the governing parameters. It is worth mentioning that these behaviors were obtained analytically and validated numerically. By focusing on the stable solution (validated numerically), it can be seen that below the critical nodes, the quantities ψ 0 , Nu, and Sh are increasing functions of Sr. Above the critical node corresponding to Sh, the latter decreases by increasing Sr. However, immediately near the critical nodes, ψ 0 and Nu exhibit a behavior similar to that of Sh and their evolutions are characterized by complex variations far above these critical nodes. Another important difference observed by comparing heat and mass transfer from Figures 6(b) and 6(c) is the fact that Nu decreases by increasing Du, while Sh increases/decreases by increasing Du for Sr < − 0.1/(Sr ≥ − 0.1). Furthermore, the most important variations observed for Nu/Sh when Sr is varied are located below/above the critical nodes. Quantitatively, at Du � − 3 for instance, by increasing Sr from − 0.9 to 0.9, Nu is multiplied by a factor of 4.1 while Sh is increased by only 11.44%. For Du � 2, by increasing Sr from − 0.9 to 0.9, Sh is divided by 3.2 while Nu undergoes a relative decrease of about 14.82%. e combined effects of Du and Sr on the temperature and concentration fields are illustrated in Figure 7 for R T � 500, N � − 0.2, Le � 1.1, Da � 0.01, Sr � ± 0.9, and Du � ± 0.9. For the combination (Du, Sr) � (0.9, 0.9), it can be seen from Figure 7 that the isotherms are very similar to the isoconcentrations. Even the quantitative comparisons (results not presented) showed that the temperature field is not very different from the concentration one. e temperature and concentration profiles obtained at midwidth of the cavity and plotted in Figure 8 show that the temperature/concentration gradients are important in the vicinity of the horizontal wall, compared to the those (i.e., gradients) around the horizontal median. By changing the combination (Du, Sr) from (0.9, 0.9) to (0.9, − 0.9), we can observe an accentuation/attenuation of the isoconcentration/isotherm distortions. e isotherms are nearly straight lines tilted with respect to the horizontal direction, which means that the temperature field in the cavity (if we except the vicinity of the short walls where the end effect is important) is quasi-linear. e temperature profile at midwidth of the cavity (Figure 8) confirms this behavior. It should be noted that the quasi-linearity exhibited by the temperature field for (Du, Sr) � (0.9, − 0.9) cannot be explained by the dominance of the diffusive regime. In fact, the diffusive regime leads to linear profiles with the isotherms being parallel to the horizontal walls. However, the isotherms obtained for (Du, Sr) � (0.9, − 0.9) are strongly tilted with respect to the horizontal direction. is means that this behavior is due to the interaction between convection, conduction, and diffusion-thermo (i.e., Dufour effect). e changes observed in the temperature/concentration field when (Du, Sr) passes from (0.9, 0.9) to (− 0.9, 0.9) are similar to those observed in the concentration/temperature field when (Du, Sr) passes from (0.9, 0.9) to (0.9, − 0.9). Hence, we can observe a quasi-linear behavior for the concentration field at (Du, Sr) � (− 0.9, 0.9). e combination (Du, Sr) � (− 0.9, − 0.9) leads to the most important variations for the temperature and the concentration among the four combinations considered in this comparison. Quantitatively, by changing (Du, Sr) from (0.9, 0.9) to (0.9, − 0.9)/(− 0.9, 0.9)/(− 0.9, − 0.9), the temperature difference between the horizontal wall at a fixed cross section passes from 0.1769 to 0.0796/0.2416/2.08. us, the temperature difference obtained with the combination (Du, Sr) � (− 0.9, − 0.9) is at least 8.6 times higher than that obtained for the other combinations. A result with the same order is also registered for the concentration difference.

Effect of Darcy Number.
e Darcy number, Da, measures the relative importance of the permeability of the porous medium as it is proportional to the latter. It is very small for the well-packed porous media and relatively large for sparsely packed ones. e effect of Da on ψ 0 , Nu, and Sh is illustrated in Figures 9(a)-9(c) for R T � 100, Le � 1.1, and different combinations (Sr, Du, N) corresponding to different regions. In these figures, it can be seen that the numerical results are in excellent agreement with the analytical ones, corresponding to the stable branches. In addition, regardless of the region considered, the parallel flow solution vanishes when Da exceeds a critical value, Da cr , that depends on the combination (Sr, Du, N). For the stable branches, the quantities ψ 0 , Nu, and Sh decrease by increasing the permeability of the porous medium. Similar behaviors were reported in previous studies [5,36]. Note that the variations versus Da of the flow intensity corresponding to the unstable branches are remarkably strong in specific ranges of the latter parameter. However, the effect of Da on Nu and Sh corresponding to the unstable branches is clearly much less important. For the combinations (0.6, − 0.5, − 2) and (0.6, 0.7, 2) corresponding to the regions 2 and 3 for which only the stable solution exists, the corresponding values of Da cr are 0.441 and 0.428, respectively. For these regions (2 and 3) Da cr can be computed analytically by solving the following equation: For the combinations (0.4, − 0.3, 2), (0.9, 0.5, 2), and (0.9, − 2, − 2), corresponding, respectively, to regions 4, 5, and 6, both stable and unstable branches exist for specific ranges of Da and the corresponding critical Darcy numbers, Da cr , marking the transition from the parallel flow towards the rest state are 0.4072, 0.1673, and 0.2938, respectively. Note that the bifurcation around Da cr is supercritical for regions 2 and 3 and subcritical for regions 4, 5, and 6 (results not reported here). For regions 4, 5, and 6, Da cr can be computed analytically by solving the following equation: A Da cr � f + 2Le 2
is means that, for Da � Da Nu�∞ , the Dufour heat flux induced by the concentration gradient compensates the ordinary heat flux induced by the temperature gradient. Such a compensation is not observed for the mass transfer as Sr ≠ Du and Le ≠ 1. For the particular combination (0.9, − 2, − 2), Sh induced by the unstable branch (in its existence range) is higher than that induced by the stable branch (see Sh behavior in the vicinity of Da cr ). Note that ψ 0 and Nu induced by the stable branches are larger than those corresponding to unstable branches for all the combinations considered in Figures 9(a)-9(b).