HYSTERESIS AND BISTABILITY BIFURCATION INDUCED BY COMBINED FLUID SHEAR-THICKENING AND DOUBLE-DIFFUSIVE CONVECTION IN SHALLOW POROUS ENCLOSURES FILLED WITH NON-NEWTONIAN POWER-LAW FLUIDS 1

,


INTRODUCTION
In recent decades, natural convection of heat and mass transfer has been characterized by fluid flows driven by buoyancy effects due to temperature and concentration gradients applied at different diffusion rates.This phenomenon is frequently encountered in nature and has many applications in the fields of science and technology, geothermal systems, heat exchangers, cooling systems, heat exchangers, microchannels, solar energy collectors, non-Newtonian chemical processes, and many others [1][2][3].Over the last few decades, most studies carried out on Newtonian double-diffusive convection in closed geometries have focused on simple geometries such as a rectangular cavity or a rectangular porous layer heated from below from the side.However, in many of the fields mentioned, the fluids used are non-Newtonian in nature.A few limited studies have been carried out on natural of power-law fluids in two-and three-dimensional enclosures [4][5][6].
Natural convection in rectangular enclosures with differentially heated vertical walls and adiabatic horizontal walls is one of the most studied configurations for Newtonian flows [7,8].In fact, Newtonian behavior at higher shear rates with dynamic viscosity dependencies as a function of temperature has been reported, as has pseudoplasticity at lower shear rates.Non-Newtonian fluids, i.e., those that go against Newton's law of viscosity, have been widely used in many industries, such as food, petrochemicals, pharmaceuticals, etc.However, when it comes to rheology, non-Newtonian fluids have been widely used.With regard to rheology, non-Newtonian fluids can be classified into Bingham plastic fluids, pseudoplastic fluids and dilatant fluids.In addition, non-Newtonian fluid models can be divided into timedependent, time-independent and viscoelastic models [9][10][11][12].Moraga et al. [13] studied the 3D natural convection heat transfer of a non-Newtonian power-law fluid in a container placed in an enclosure.The container was filled with a non-Newtonian power-law fluid and surrounded by air.Their results included isotherms and streamlines at different values of Ra for each Newtonian and non-Newtonian power-law fluid.Jahanbakhshi et al. [14] studied the influence of the magnetic field on the natural convection of a non-Newtonian fluid in an L-shaped cavity.They found that the rate of heat transfer increased for shear-thinning fluids and decreased for shear-thickening fluids.They also found that the Nusselt number decreases with increasing Pr when the fluid shear-thins.Lamsaadi et al. [15] have numerically simulated the natural convection of a non-Newtonian fluid in a vertical cavity.Their results show that a shear-thinning fluid contributes to fluid flow and a shear-thickening fluid reduces it.Shear-thinning fluid also increases heat transfer.Makayssi et al. [16] presented an analytical and numerical experimental study of natural double-diffusion convection in a shallow twodimensional horizontal enclosure.These authors concluded that the characteristics of fluid flow, heat and mass transfer appear to be highly sensitive to the flow behavior index, n.Thus, compared to the Newtonian case (n = 1), shear-thinning
behavior (0 < n < 1) enhances fluid flow as well as convective heat and mass transfer.Kefayati [17] considered non-Newtonian fluid in a square enclosure and discussed the effect of Reynold's and magnetic field inclination angle on fluid flow.Zyła et al. [18] carried out a rheological study on nanofluids containing a mixture of nanodiamond and graphite with different ash fractions.They also studied a mixture of nanodiamond and graphite with different ash fractions.Non-Newtonian behavior and viscoelastic structure were also observed.The transient buoyant convection of a non-Newtonian fluid in a square enclosure, heated differentially from the vertical, has been studied numerically by Kim et al. [19].Getachew et al. [20] have described the flow type and heat and mass transfer patterns of double-diffusive inside a porous enclosure saturated with a non-Newtonian fluid of power using scaling arguments, verifying the results obtained by experimental.Moradi et al. [21] have carried out an experimental study of heat transfer rate enhancement of multi-walled carbon nanotube nanofluids in a double-pipe heat exchanger with aluminium porous media.Jena et al. [22] have studied the double-diffusive free convection of a non-Newtonian Ostwald-De Waele fluid inside a 2D cavity with partially active vertical walls.Lounis et al. [23] studied the impact of Dufour and Soret effects on double-diffusive convection in an inclined square enclosure using the Carreau-Yasuda model to model the rheological behavior of a non-Newtonian fluid.
Their main conclusion is that for different values of the power-law index, the Lewis number increases heat and mass exchange.Rebhi et al. [24,25] investigated the convection generated in porous cavities containing a binary fluid and discovered that the drag parameters had a significant impact on the stability and onset of subcritical and Hopf bifurcation.Rebhi et al. [26] studied Rayleigh-Bénard thermosolutal convection instabilities in shallow enclosures with finite aspect ratios that influenced by the rheological behavior of non-Newtonian fluids.They showed that the bistable convection regime exists for every aspect ratio of the enclosure, regardless of the type of thermal and solutal boundary conditions.Bihiche et al. [27] analyzed the effect of buoyancy ratio on natural convection and on heat and mass transfer rates in a horizontal rectangular enclosure filled with a non-Newtonian power-law fluid.The emergence and development of double-diffusive convection for both aiding and opposing flows were also studied.Consequently, natural convection in rectangular/square cavities filled with porous layers saturated with a non-Newtonian fluid has been studied by many researchers [28][29][30].Nield et al. [31] studied the effect of high heterogeneity on the occurrence of vertical density gradient-induced convection in soil.On the occurrence of convection induced by a vertical density gradient in a saturated porous medium governed by Darcy's law.Devakar et al. [32] have simulated the flow of non-Newtonian fluids in a straight, uniform square duct through a porous medium.They observed that the velocity and volumetric flow rate decrease with increasing torque stress and porosity parameters.Uphill et al. [33] studied the flow of nanofluids in porous media.They showed that nanofluids containing particles smaller than the size of water are more efficient in porous media.Lakshmi et al. [34] studied the effect of Soret and Dufour diffusion on natural convection in a saturated porous medium.The results indicate that the Nusselt number increases linearly with increasing Dufour parameter to facilitate buoyancy.Numerical study of Soret and Dufour effects on the natural convection flow of a vertical plate with power-law heat flux embedded in a porous medium is studied by Tsai and Huang [35].Zhu et al. [36] have numerically studied the natural double-diffusive convection of power-law fluids in a porous cube using a generalized non-Darcy model.They found that the impact of different power law indices on convection is mainly manifested in rheological properties, showing that shear-thinning fluids are more efficient at enhancing heat and mass transfer than shear-thickening fluids.Madhua et al. [37] investigated new entropy generation characteristics in the fully developed heat transport of a non-Newtonian Carreau fluid in an inclined microchannel.They found that the velocity is maximal at the center of the microchannel when using increasing values of Grashof number Gr, while entropy production is maximal at the center of the channel.Ben Khelifa et al. [38] have derived an analytical solution to characterize the onset of motion, heat and mass convection for a binary power law fluid in a shallow porous enclosure.The occurrence of nonlinear convection in a horizontal porous layer saturated with a shear-acting fluid is studied by Bensilakhal et al. [39].The most interesting results of their study are the demonstration of the existence of a bi-stability phenomenon, i.e., the existence of two steady-state solutions, which had not previously been observed in the convection of non-Newtonian fluids.
Many studies have focused on the effect of the shear-thinning of non-Newtonian fluids, which attracts more attention when it comes to heat transfer problems.The Marangoni in a shallow rectangular cavity of a power-law fluid was studied by Naimi [40].Chen [41] examined the influence of Marangoni convection for a power-law liquid film on an unsteadily stretching aluminium foil.Lamsaadi [42] studied the natural convection of non-Newtonian power-law fluids subjected to transverse temperature gradients inside a horizontal rectangular space.These authors noted the combined effect of heat flow ratio, power-law index and Rayleigh number on flow intensity and heat transfer rates.Lamsaadi et al. [43] have shown that the depends only on the nominal Rayleigh number Ra and the power-law index, n, for large values of the aspect ratio and the nominal Prandtl number Pr. Habibi et al. [44] described the natural laminar convection of a non-Newtonian fluid between two horizontal eccentric square ducts under equilibrium conditions.Recently, Alloui et al. [45] analyzed the onset of fluid motion for power-law fluids.It was shown that, for shear-thinning fluids, the onset of convection is subcritical.For shearthickening fluids, convection occurs at a supercritical Rayleigh number of zero.Khali et al. [46] studied double-diffusive convection in a power-law fluid.The results indicate that the fluid structure is more important for thermal base flow.Ohta et al. [47] have numerically investigated transient heat transfer in a square cavity heated from below and cooled from above.For shear-thinning fluids, the Sutherby model was used.Their study also reveals an increasing magnitude of shear-thinning, based on a transient analysis of the natural convection of shear-thinning-fluids in rectangular enclosures with differential heating.Hojjat [48] described Reynolds and Prandtl numbers for non-Newtonian fluids; the coefficient index and power-Hysteresis and Bistability Bifurcation Induced by Combined Fluid Shear Thickening...

EEJP. 1 (2024)
law index are very important in determining the Reynolds number and initial velocity of the fluid at the tube inlet.Reynolds and Prandtl numbers are described as follows for non-Newtonian fluids.Solomatov and Barr [49,50] have numerically investigated the occurrence of convection in non-Newtonian power-law fluids under the effect of temperature-dependent viscosity.The authors summarize their results in terms of simple scaling relations.The primary focus of this investigation centered on the phenomenon of thermosalutal convection within a horizontal porous layer filled with a non-Newtonian fluid, subjected to heating and salting from below.The non-Newtonian fluid behavior was accurately characterized through the implementation of the power-law model, accepting a diverse range of fluid types, including shear-thinning ( < 1), shear-thickening ( > 1), and Newtonian fluids ( = 1).Employing a time-accurate finite difference method, the comprehensive nonlinear governing equations were numerically resolved.Additionally, an analytical solution for shallow enclosures was derived.The outcomes of this study significantly contributed to an enhanced comprehension of the impact of various governing parameters on diverse convective bifurcations.Notably, the shapes of bifurcation branches underwent substantial alterations, particularly in scenarios involving opposing double-diffusive flows.The most noteworthy revelation was the emergence of bistability bifurcation with increasing shear-thickening behavior of the fluid, leading to the occurrence of two distinct and stable solutions under identical flow conditions, a phenomenon hitherto unobserved in the context of Newtonian fluids.

PROBLEM DEFINITION AND MATHEMATICAL FORMULATION
The problem under study is sketched in Figure 1.It's a tow-dimensional horizontal shallow cavity filled with non-Newtonian binary fluids saturated porous media.Choose an appropriate coordinate system where the  ,  are the longitudinal and transversal axis, respectively.The enclosure is of height  and length  .The short vertical end walls are thermally insulated, while the horizontal walls are submitted to constant heat flux,  , and solutal flux,  .In accordance with previous statements made by Pascal [51,52], the model of laminar flow of a non-Newtonian power-law fluid through a porous medium is explained as follows: ( ) , and 8 1 3 (2 ) where  is the superficial velocity,  and  the porosity and the permeability of the porous medium respectively,  , the apparent viscosity, ℎ the consistency index and  the power-law index.In the above model the rheological parameters ℎ and  are assumed to be temperature independent.
The equations governing the heat flux  and matter flux  , with respect to the thermal and solute gradients in the binary fluid mixture, are expressed by De Groot and Mazur [53], as follow: Herein, where ′ and ′ represent respectively the temperature and the concentration of the fluid at a point in the system,  and  represent the thermal conductivity and molecular diffusion coefficient of the species, respectively.
If we make the assumption that the flowing fluid and the porous medium are in local thermodynamic equilibrium throughout, with constant properties for both the fluid and the porous medium, and apply the Boussinesq approximation, which has been used in the past by many authors Amari et al. [29], Bian et al. [54,55], and Benhadji and Vasseur [28].The equations that govern the current problem can be expressed as follows: . 0 where  the density of the non-Newtonian fluid. ,  and  represent the temperature, concentration, and mass density of reference,  the gravitational acceleration, () and () are the heat capacities of fluid and saturated porous medium, and  = /() the thermal diffusion coefficient. ,  are respectively the thermal and solutal expansion coefficients.They are defined by , , 0 0 1 1 , By using  , / ,  =   /, and  =  /Δ′ to scale length, velocity, temperature, and concentration, respectively, and eliminating the pressure term from Equation ( 5) through standard methods, it can be demonstrated that the dimensionless governing equations can be formulated as follows The dimensionless parameters who were extracted from Eqs. ( 9)-( 11) are  = ′  / is a thermal Rayleigh number, the Lewis number  =   ⁄ , the buoyancy ratio ⁄ is the heat capacity ratio, and Ψ is a dimensionless stream function defined as , u v dy dx Remember that it is possible to improve the Newtonian expressions by simply set =1 and  = 1.The nondimensional boundary conditions at the enclosure walls are as follows where  =  / is the cavity aspect ratio.The problem definition is rendered complete by the simultaneous consideration of Eqs. ( 9)-( 13) along with the boundary conditions, as denoted by Eqs. ( 14) and (15).It is important to note that the solution to this problem is intrinsically dependent on the values of the parameters  , , ,  and .
From an engineering perspective, the focus lies in the calculation of heat and mass transfer rates expressed through local and average Nusselt and Sherwood numbers, denoted as  ,  , ℎ , and ℎ , respectively.In the current notation, the computation of  ,  , ℎ , and ℎ is performed as follows: where the integral is performed using Simpson's method.

NUMERICAL SOLUTION
For simple geometry problems Fig. 1, the finite difference method used to numerically solve the governing equations ( 9)- (11).The energy and concentration equations are discretized with a second-order centered scheme.For each time step the implicit method with alternating directions (ADI) gives rise to two tri-diagonal matrix systems to be solved, one Hysteresis and Bistability Bifurcation Induced by Combined Fluid Shear Thickening...

EEJP. 1 (2024)
resulting from the implicit discretization in  and the other of the implicit discretization in .The (ADI) method divides the time step into two equal parts, in the first half-time step, the system is implicitly in  but explicitly in  and in the second half-time step, the system is implicitly in  and explicitly in .For each time step, the solution is obtained by scanning the computational domain in the  direction then in -direction.Knowing the fields of temperature and concentration at time ( 1), to solve the discretized motion equation at each time step successive over relaxation method (SOR) is used which is an explicit method directly giving the value of  to the instant ( 1).The convergence criterion for solving Eq. ( 9) is provided.
, 10 The selection of mesh size is a careful consideration aimed at striking a suitable balance between computational efficiency and result accuracy.The approach involves grid refinement iteratively until the numerical solution converges, with a reasonable accuracy threshold, to the analytical solution developed in the subsequent section.This convergence is achieved under the prescribed values for , , , , and  .As indicated in Table 1, a uniform grid of dimensions 221x141 proves to be well-suited for accurately modeling the flow, temperature, and concentration fields within a cavity between 6 and 10.The chosen time-step sizes for simulations range between 10 -6 and 10 -3 .apparent viscosity lines are uniformly distributed with specified increments  between the minimum value,  , on the center of cavity and the extremum value,  , on the center when  = 0.6, this is the opposite of  = 1.4,and it's clearly constant ( = 1) when  = 1.The numerical findings depicted in Figure 2 validate the theoretical assumption made in the analytical solution, indicating that the flow patterns within the central region of the cavity remain parallel regardless of the power-law index, .It is important to note.However, that the scope of this current investigation is confined to the examination of unicellular convection exclusively.
To validate the accuracy of the current numerical solutions, Tables 2 and 3 presents the obtained numerical results, including the center stream function,  , value, as well as the local Nusselt, , and Sherwood, ℎ, numbers.A comprehensive comparison with the findings reported by Amari et al. [29] and Ben Khelifa et al. [38] reveals a very good level of agreement.

ANALYTICAL SOLUTION
An analytical solution for the full governing equations, ( 9)-( 11) is in general impossible, except for certain cases and assumptions for which the equations simplify considerably.In the case of a porous layer of large aspect ratio ( ≫ 1) with active walls exposed to constants fluxes of heat and mass, it is possible to find an approximate analytical solution.This last, obtained using the concept of parallel flow, the generated flow becomes relatively parallel along the long walls of the cavity, which of course makes it possible to neglect the component of the velocity perpendicular to these walls.
The streamlines, at the center of the cavity, become parallel to the x-axis.In other words, the current function Ψ becomes a function of the y ordinate only.We can then write: The temperature and concentration profiles are given by the sum of two terms, the first one defining a linear longitudinal variation and the other giving the transverse distribution: where  and  are constants expressing temperature and concentration gradients along the x-direction.Substituting approximations (18) and (19) in the governing equations ( 9)-( 11) and assuming steady state of flows, the heat and mass fluxes, we get the following ordinary differential equations: Solution of Eqs. ( 20)-( 22) satisfying the boundary conditions given by Eq. ( 14) are: ( ) [ ] ( ) ( ) [ ] ( ) ( ) [ ] ( ) where: The expression of the velocity distribution  along y derived from the stream function is: The concept of parallel flow loses its viability due to the turning flow at the end of perpendicular walls.However, in the case of thermal and mass conditions in the end walls, the value of constant temperature and concentration gradients,  and  respectively are determined by considering the arbitrary control volume of Fig. 1, the energy and material balances in the control volume at any -direction are written according to the following forms: By substituting the temperature, concentration and velocity Eqs. ( 24)-( 26) into Eqs.( 27), we get the following expression: Equations ( 28) and ( 29) has been solved numerically for  and  using the method of Newton-Raphson.From Eqs. ( 24), ( 25) and ( 16), the Nusselt and Sherwood numbers are given by: [ ] ( ) [ ] ( )

LINEAIR STABILITY ANALYSIS
In this section, an investigation into the stability of the steady state is conducted.The overall convective solution encompasses a foundational solution (Ψ ,  ,  ) with a perturbation solution (Ψ ,  ,  ).
In the context of an infinite horizontal fluid layer, the perturbation solution for the stream function, temperature, and concentration field can be formulated as follows:

, pt i x p pt i x p pt i x p x y t e F x y T x y t e G x y S x y t e H x y
In this formulation, (, ), (, ), and (, ) represent the spatial profiles of the perturbation for the stream function, temperature, and concentration, respectively. characterizes the growth rate of the perturbation,  is a real number which designates the wave number.The variables Ψ ,  , and  denote infinitesimal unknown perturbation amplitudes for the stream function, temperature, and concentration.
In order to assess the stability of the rest state solution, a tiny perturbation is introduced, allowing the complete solution to be expressed as follows: EEJP. 1 (2024) Saleh Khir, et al.

p b p x y t x y x y t T x y t T x y T x y t S x y t S x y S x y t
The apparent viscosity, denoted as  , can be further decomposed as: The linearized form of the perturbed apparent viscosity is expressed as follows: ( ) Upon substituting equations ( 32)-( 35) into equations ( 9)-( 11) and neglecting the second-order nonlinear terms, the resulting linearized stability equation is as follows:

(
) ( ) where: By employing the variational formulation through the Galerkin technique, the integrated linear equations are derived as follows: The constants  Hysteresis and Bistability Bifurcation Induced by Combined Fluid Shear Thickening...

RESULTS AND DISCUSSION
The results of this study focus on the effect of the power-law index, , the Rayleigh number,  , the buoyancy ratio, , and the Lewis number, , on the bistability convective flows, heat and mass transfer behaviors, and on the thresholds of the onset of various convective bifurcations in a porous layer, as predicted by the analytical solution.The findings have been confirmed through numerical analysis of the complete governing equations.The primary revelation in this research involves identifying a convective bistability phenomenon occurs from the interactional effect between the slower diffusing component and the shear-thickening fluids ( > 1) properties.The study encompasses both numerical and analytical outcomes, which fall within the following typical ranges: 0.6 to 1.4 for ',' 1.0 to 20.0 for ',' and -0.2 to 0.2 for '.' Our research is primarily concerned with the onset of bifurcation.Consequently, we restrict our analysis to Rayleigh numbers  of relatively modest magnitudes, specifically  ≤ 100, in close proximity to the critical points.Furthermore, in order to establish the presence of bistability in convection across various aspect ratio and thermal and solutal boundary conditions, we systematically vary the parameter '' within the range of 6 to 10, while adopting Neumann boundary conditions.This approach serves to generalize the conditions under which bistability in convection phenomena can manifest.Our initial investigation focuses on a shallow enclosure, benefiting from the availability of an analytical solution to provide guidance in this context.
Figures 3(a   Notably, for a Newtonian fluid, the velocity distribution follows a linear pattern, in agreement with the findings of Amari et al. [29] and Vasseur et al. [56].The depiction of apparent viscosity in Figure 3(b) elucidates a diminishing trend as the power-law index, , decreases.This signifies a reduction in the significance of convective motion when the fluid exhibits shear-thickening behavior ( > 1).Conversely, contrasting effects are observed for shear-thinning fluids ( < 1), as evidenced in Figure 3(b).These observed phenomena align with the findings reported by Amari et al. [29], Dharmadhikari and Kale [57], as well as Chen and Chen [58,59], and it's clearly constant ( = 1) in the case of Newtonian fluid ( = 1).In Figures 3(c) and 3(d), it is observed that the temperature and concentration profiles exhibit an augmentation with increasing values of the power-law index, .Additionally, all curves maintain a constant slope at  = −0.5 and  = 0.5, attributed to the imposition of a constant heat and mass flux on the horizontal walls.
As an introductory exploration into the identification of bistability in convection, we present a graphical representation with an algorithmic scale in zoomed image Figure 4, illustrating the bistability bifurcation concerning flow intensity as a function of the Rayleigh number.The figure encompasses three distinct Rayleigh number thresholds, each bearing significance.Specifically, it is well-established that, in the context of a negative buoyancy ratio ( < 0), a threshold similar to that of subcritical flows emerges, denoted as  .This threshold may manifest either below or above the supercritical threshold,  .Furthermore, an additional threshold, referred to as  , consistently surpassing  , delineates a hypothetical backward bifurcation, symbolized by a dashed-dotted arrow.This backward bifurcation exhibits a transition from low to high finite amplitude convective states, forming an "S"-shaped bifurcation pattern.
The phenomenon of bistability bifurcation, achieved through systematic variations in Rayleigh numbers in both ascending and descending sequences, engenders a hysteresis loop.Within this loop, high and low convective states persist concurrently under identical flow boundary conditions and governing parameters.
In Figure 4, we depict three distinctive bistability regimes.The first of these regimes, referred to as the subcritical regime, manifests when the subcritical threshold  is significantly lower than the supercritical threshold  , specifically when, , is less than 1.107.The second regime, known as the transcritical regime, materializes when  equal  , denoting a condition where, , equals 1.107.The third regime, termed the supercritical regime, emerges when  considerably exceeds  , signifying, , values greater than 1.107.Across all three of these regimes, the characteristic "S" shaped bistability bifurcation remains intact, resulting in the coexistence of two distinct and stable convective states.Notably, the bistability in convection ceases to exist when the lower and upper bifurcation points,  and  , converge (i.e., when  equal 1.55), giving rise to a well-defined inflection convection state characterized by a sharp increase in heat and solute transfer rates.
Table 4 provides a comprehensive presentation of the critical Rayleigh numbers denoting the onset of motion ( ,  and  ) as a function of the power-law index, .A discernible trend emerges from the table, underscoring the notable influence of the power-law index, particularly in the context of thickening fluids.It is remarkable that such an influence does not manifest in the case of shear-thinning fluids, highlighting a significant distinction in the behavior of bistability convection.
In Figure 5, characteristic bifurcation diagrams are presented, illustrating the dependencies of Ψ , , and ℎ on  and  under the specified conditions of  = 5 and  = −0.1.The outcomes encompass a range of power-law index, spanning from  = 0.6 to 1.4.The curves featured in the plots represent the predictions derived from the present analytical/numerical nonlinear models.The solid lines denote the stable branches, while the long-dashed lines signify unstable branches.The numerical solutions derived from the full governing equations are represented by solid circles.Notably, Good agreement is observed between the results of these two nonlinear theories.The profile obtained for a Newtonian fluid  = 1.0 corresponds to convection induced by the imposition of a constant heat and mass flux on the horizontal boundaries of the system.This configuration leads to the emergence of a saddle-node bifurcation at a subcritical Rayleigh number, denoted as  =15.63.The determination of the subcritical Rayleigh number involved a numerical search within the analytical model, where the value of  was identified such that the inverse of the derivative of  with respect to  equaled zero.With an increasing in the power-law index, denoted as , the bifurcation behavior undergoes a notable transformation contingent upon the magnitude of power-law index.Across the range of cases spanning from 1 <  < 1.3, the bifurcation curve exhibits the distinctive presence of two turning saddle-node points ( and  ).These points connect the two stable branches, delineated by an intervening unstable branch represented by long-dashed lines.The upper stable branch, originating from  , aligns with the conventional subcritical bifurcation observed earlier for  = 1.Conversely, the lower stable branch, commencing at  , corresponds to a supercritical bifurcation.However, it is noted that the latter exist within the interval of  ≤  ≤  .Under these circumstances, a bistability regime manifests EEJP. 1 (2024) Saleh Khir, et al.
itself between the two thresholds ( ≤  ≤  ).Beyond the supercritical Rayleigh number,  , the outcomes derived from the numerical solution of the full governing equations suggest that, when initiating computations from the initial rest state conditions, the solution traces a hysteresis loop as illustrated by arrows in the zoom included in Figure 5(a.2;b.2 and c.2). Corresponding heat and solute transfer rates are presented in Figures 5(b) and 5(c).Under the specified conditions of  = −0.1 and  = 5, the presence of bistability bifurcation is confirmed within the power-law index range of 1 <  < 1.3.Notably, for  ≤ 1, the observed bifurcation assumes a subcritical.Within the zone of bistability, Figure 6 showcases the solutions at  = 18, obtained through the numerical resolution of the comprehensive governing equations.These visual representations encompass streamlines, isotherms, isoconcentrates, and apparent viscosities, presented in a sequential order from top to bottom.The solution portrayed in Figure 6(b), attributed to the lower branch, originates from the initiation of numerical calculations using the initial rest state conditions.However, the solution related to the upper branch, Figure 6(a), is exclusively attained by following the delineated hysteresis path on the curve for  = 1.1, with an increase in the power-law index, the region of bistability contracts, as evidenced by the outcomes derived for  = 1.1, and diminishes entirely as the fluid became more and more shear-thickening, exemplified in instances with  = 1.3.In such circumstances, a resultant pitchfork-type bifurcation emerges, characterized by the transition from the rest state to a convective regime, occurring at a supercritical Rayleigh number,  .
Table 5 gives the critical values of the Rayleigh numbers characterizing the onset of motion ( ,  , and  ) as a function of the power-law index, .This table clearly indicates that bistability convection was significantly affected by the power-law index, .In Figure 7, a stability diagram is presented to illustrate the influence of the power-law index denoted as, , and the buoyancy ratio, , on critical Rayleigh numbers that govern the onset of subcritical  and supercritical  convection, as well as an additional critical point denoted as  , which corresponds to a turning saddle-node point where a backward bifurcation occurs, resulting in a transition from lower to higher convective states.These findings are presented within the context of  = 5.0.Upon initial examination, it becomes evident that  ,  and  exhibit significant decreases as the buoyancy ratio, , increases.This decrease can be attributed to the collaborative influence of shear-thickening behavior and the slower diffusing solute, which collectively augment the subcritical convection phenomenon.A more comprehensive depiction of the impact of buoyancy ratio, , and Rayleigh number,  , on  , , and ℎ is presented in Figure 8 for  = 5 and  = 1.2.The outcomes are derived for varying buoyancy ratio values from 0.2 to -0.2.As evident from Figures 8(a)-8(c), for a given Rayleigh number, both the intensity of convection,  , and the ensuing heat, , and solute, ℎ, transfer rates exhibit a decline as the buoyancy ratio, , diminishes.The plot notably illustrates the influence of  on the existence of various bifurcation curves, as elucidated in Figure 8. Consequently, a decrease in the buoyancy ratio, , leads to the evolution of bifurcation curves from subcritical to supercritical behavior.Overall, the influence of  on the onset of motion ( ,  , and  ) is elucidated in Figure 7.Under the specified conditions of a power-law index,  = 1.2, and a Lewis number,  = 5, the presence of bistability bifurcation is confirmed within the buoyancy ratio range of  =-0.18 to -0.082, (i.e., under opposing flow conditions  < 0), Notably, for  = −0.2, the observed bifurcation assumes a subcritical.In the region characterized by bistability, the solutions corresponding to  = 22, derived from the numerical solution of the comprehensive governing equations, are depicted in Figure 9.These graphical representations sequentially present streamlines, isotherms, isoconcentrates, and apparent viscosities from top to bottom.The solution illustrated in Figure 9(b), representing the lower branch, is obtained by initializing the numerical computations with the rest state as the initial conditions.Conversely, the solution for the upper branch, as shown in Figure 9(a), is exclusively acquired by tracing the hysteresis path delineated on the curve for  = −0.14.As the buoyancy ratio becomes large, the bistability region shrinks, as shown by the results obtained for  = −0.14, and it completely disappears when the buoyancy ratio becomes very large, as displayed for thermal convection  = 0 and for the aiding flow ( > 0).For this situation, the resulting pitchfork type bifurcation, characterized by a transition from the rest state to a convective regime, takes place at a supercritical Rayleigh number,  .Table 6 provides the critical Rayleigh numbers, delineating the onset of motion ( ,  , and  ), as functions of the buoyancy ratio, .This table clearly indicates the impact of the buoyancy ratio, , on bistability convection.In Figure 10, a stability diagram is presented to illustrate the influence of the power-law index denoted as, , and the Lewis number, , on critical Rayleigh numbers discus in figure 7.These findings are presented within the context of  = −0.1.Upon initial examination, it becomes evident that the bistability region presented in zoomed image Figure  number progressively increases, the bistability region manifests and undergoes significant expansion, persisting within the range below 10, indicative of thermal diffusion dominant over mass diffusion.The graphical representation in Figure 11 provides insight into the notable impact of the magnitude of the Lewis number on the characteristic "S"bifurcation shape within the bistability region.Specifically, the "S"-bifurcation curve becomes increasingly conspicuous with higher Lewis numbers (i.e.,  = 7 and 9).At the subcritical bifurcation point, the flow exhibits finite amplitude, marking a discernible transition from a rest state to a convective state.This transition is characterized by a sudden and substantial augmentation in flow intensity, Nusselt, and Sherwood numbers.Such enhancements are attributed to a pronounced reduction in the apparent viscosity, as depicted in Figures 11(a)-11(c).Table 7 shows the impact of the Lewis number, , on the critical Rayleigh numbers ( ,  , and  ) for  = 1.1, and  = −0.1.It was found that the bistability regime occurs when Lewis number, , values are within the domain ]3.4;10[.

CONCLUSIONS
In this paper, the problem of thermosalutal convection of power-law fluid saturated porous media contained in a horizontal enclosure subject to vertical constant fluxes of heat and solute, has been investigated analytically and numerically.The effect of the thermal Rayleigh number  , the power-law index , buoyancy ratio , Lewis number , and the aspect ratio of the cavity, , on the onset of linear and nonlinear convective motion, as well as heat and mass transfer rates.for the case of a shallow enclosure  ≫ 1, an analytical solution, based on the parallel flow approximation, has been formulated.The results show good agreement between the numerical and analytical solutions, and a strong influence of the governing parameters on the onset of convective motion and the resulting convective heat and mass transfer rates.
In the first part, results are presented for the velocity, apparent viscosity, temperature and concentration fields, revealing a significant sensitivity to the power law-index n.It is shown that an increase in the power-law index (dilatant fluids) enhances the apparent viscosity while a decrease (pseudoplastic fluids) reduces it, and it's constant for the case of Newtonian fluid ( = 1), it is observed that the temperature and concentration profiles exhibit an augmentation with increasing values of the power-law index, .
It was observed that the intensity of convection Ψ and the resulting heat  and mass ℎ transfer rates experienced enhancement with a decrease in the value of , and , while reduces it when  decrease.
The resulting nonlinear solution indicated that the threshold of finite amplitude motion occurred at a supercritical Rayleigh number ( ), which considerably decreased upon increasing the power-law index , and , and  decrease.
Specifically, the manifestation of bistability convection phenomena, denoted by the coexistence of two distinct stable solutions under identical flow and boundary conditions, is contingent upon specific values of , , and .Notably, this phenomenon is associated with the characteristics of a dilatant fluid ( > 1) and opposing buoyancy forces ( < 0).The bistability bifurcation was observed to traverse three saddle-node points, corresponding to the thresholds  ,  , and  .This bifurcation exhibited three branches, resulting in an "S" shape curve.The first branch bifurcated forward from  to  , the second bifurcated backward from  to  but was deemed unstable, and a third stable branch bifurcated forward from  .Furthermore, it was demonstrated that this bistability phenomenon could manifest irrespective of the aspect ratio of the enclosure and the types of thermal and solutal boundary conditions.
In the case of opposing convective flows, there was a collaborative influence of both shear-thinning and slowdiffusing solute effects, resulting in the augmentation of subcritical convective flows.Conversely, for aiding convective flows, a cooperative impact of shear-thickening and slow-diffusing solute effects was observed, leading to the enhancement of supercritical convective flows.
The occurrence of bistability convection was attributed to the cooperative counteracting influences within the fluid system.Specifically, the dilatant fluid.Simultaneously, the slower diffusing component in the binary mixture tended to maintain uniformity within the core region of the system while establishing a pronounced solute gradient near the walls.

Figure 1 .
Figure 1.Sketch of the enclosure and coordinates system.

Figure 2 .
Figure 2. Contours of streamlines, , Temperature, , Concentration, , horizontal velocity, , and apparent viscosity,  , for  = 100,  = 5,  = −0.1 and various values of  Figure 2 present the contours of streamlines, , temperature, , concentration, , horizontal velocity, , and apparent viscosity,  , obtained numerically, for two different liquids with power-law index of n = 0.6 and 1.4 and for a Newtonian fluid (n = 1), at  = 100,  = 5, and  = −0.1.The streamlines are uniformly distributed with designated increments, , spanning from zero at the boundaries to the maximum value,  , at the center.The isotherms and isoconcentration lines are uniformly distributed between the successive positions of the peak temperature  and maximum concentration  situated lower edge of the right vertical wall.On the other hand, the minimum temperature and minimum concentration  are situated along the upper edge of the left vertical wall.Notably, the dimensionless temperature and concentration at the cavity's center are both normalized to zero.Furthermore, owing to the centrosymmetry inherent in the problem, the minimum and maximum temperatures, as well as concentrations, exhibit equivalent magnitudes but with opposite signs.The horizontal velocity lines are uniformly distributed from the maximum velocity,  , located on the lower horizontal wall and the minimum velocity,  , on the upper horizontal wall.The
and basic apparent viscosity, respectively.
)-(d) show the horizontal velocity, , apparent viscosity,  , temperature, , and (d) concentration, , profiles at the center of the layer ( = 0) for  = 100,  = 5,  = −0.1, and various values of . the results are presented for −0.5 ≤  ≤ 0.5.The analytical solution, represented by solid lines, aligns closely with the numerical results, denoted by solid circles, demonstrating a favorable agreement between the present analytical solution and the numerical solution.

Figure 3 (
a) reveals a discernible decrease in velocity with increasing values of the power-law index, , and highlights that the maximum velocity values occur at the boundaries.This behavior stems from the modeling of the porous medium in accordance with Darcy's law, permitting fluid slip at the solid boundary.

Figure 4 .Figure 5 .
Figure 4. Bifurcation diagram in terms of Ψ as the function of  and
decrease as the Lewis number, , decrease until a value of  = 2.4, where  = = and  = 1, when the values of Le exceeds approximately 10 the bistability regime disappear the invers of Rebhi et al[26].For this interval of Lewis number, , and with an increasing of , the domain of  expands and remains above1.

Table 6 .
Dependence of  ,  and  on  for  = 5.0 and  = 1.2.