Natural Convection Flow over a Vertical Permeable Circular Cone with Uniform Surface Heat Flux in Temperature-Dependent Viscosity with Three-Fold Solutions within the Boundary Layer

The aim of this study is to investigate the effects of temperature-dependent viscosity on the natural convection flow from a vertical permeable circular cone with uniform heat flux. As part of numerical computation, the governing boundary layer equations are transformed into a non-dimensional form. The resulting nonlinear system of partial differential equations is then reduced to local non-similarity equations which are solved computationally by three different solution methodologies, namely, (i) perturbation solution for small transpiration parameter (ξ), (ii) asymptotic solution for large ξ, and (iii) the implicit finite difference method together with a Keller box scheme for all ξ. The numerical results of the velocity and viscosity profiles of the fluid are displayed graphically with heat transfer characteristics. The shearing stress in terms of the local skin-friction coefficient and the rate of heat transfer in terms of the local Nusselt number (Nu) are given in tabular form for the viscosity parameter (ε) and the Prandtl number (Pr). The viscosity is a linear function of temperature which is valid for small Prandtl numbers (Pr). The three-fold solutions were compared as part of the validations with various ranges of Pr numbers. Overall, good agreements were established. The major finding of the research provides a better demonstration of how temperature-dependent viscosity affects the natural convective flow. It was found that increasing Pr, ξ, and ε decrease the local skin-friction coefficient, but ξ has more influence on increasing the rate of heat transfer, as the effect of ε was erratic at small and large ξ. Furthermore, at the variable Pr, a large ξ increased the local maxima of viscosity at large extents, particularly at low Pr, but the effect on temperature distribution was found to be less significant under the same condition. However, at variable ε and fixed Pr, the temperature distribution was observed to be more influenced by ε at small ξ, whereas large ξ dominated this scheme significantly regardless of the variation in ε. The validations through three-fold solutions act as evidence of the accuracy and versatility of the current approach.


Introduction
Natural convection occurs in the environment. Different industrial applications, closed containers, or any type of chamber, such as a greenhouse, are also good examples of locations where convective heat transfer occurs. In order to characterize different heat transfer applications, heat flux is one of the major indices to be considered. Uniform heat flux is generation of fluids within the incompressible boundary layer to understand the changes in velocity and concentration profiles, which could only be explored due to the consideration of temperature-dependent viscosity. Gladys and Reddy [32] also accentuated the role of the temperature-dependent viscosity of non-Newtonian nanofluids through the accelerating vertical plate to analyze the nonlinear buoyancy impacts.
In the literature, most of the works consider the constant viscosity of the fluid within the boundary layer and failed to prove the accuracy of the model in various ranges of Pr numbers. The current study considers both small (for example, 0.05) and large (for example, 0.7) Pr numbers. Therefore, it is essential to validate the current approach with different numerical parametric tests to gain more confidence in the accuracy and versatility of the approach. Therefore, more than one numerical solution/validation should be conducted to gain a better understanding of the accuracy. Among different numerical techniques, the Keller box method has been a popular and widely accepted numerical technique for nearly five decades [33], with more applications being added to this scheme recently. For example, Kamran et al. [34] have considered the Keller box approach to describe the Jefferey-Hamel flow by considering different non-dimensional parameters to obtain solutions. Reddy et al. [35] also obtained implicit finite difference results by the Keller box technique to study the Joule heating and associated chemical reactions on the magneto Casson nanofluid. Therefore, another evidence of expanding applications of the Keller box method was duly noted. On the other hand, perturbation and asymptotic solutions were also found to be highly efficient for small and large transpiration parameters (ξ), respectively, for the last few years [36]. Therefore, the validations through a combination of conventional and modern approaches should be adequate to prove the accuracy of the current approach and will provide a new dimension to future study.
The present study aims to investigate natural convection flow in a vertical permeable cone with uniform heat surface flux as viscosity varies in the computational model. The present model is also valid for fluids with a low Prandtl number (Pr), thus resolving the shortcomings of most of the literature in the past decades. The non-dimensional viscosityvariation parameter (ε) has been included in the model, along with the pseudo-similarity variable (η), to tune the fluid characteristics in the model. The suction parameter (ξ) has been included at various ranges to observe the effect of changes in the rate of heat transfer and shear stress. The sensitivity analyses have been conducted with both constant and variable Pr numbers to observe the behaviors of the temperature and velocity distribution of the fluid within the boundary layer. The model has been validated with both fixed and variable Pr numbers to showcase the uniformity as well as the accuracy of the approach. The three-fold numerical solutions have been conducted by implementing the aforementioned parameters. The implicit finite difference method together with the Keller box scheme for all ξ has been conducted, followed by perturbation for small ξ and, finally, the asymptotic solution for large ξ. To the authors' knowledge, there has not been any published work that considers all the characteristics and scopes mentioned above.

Formulation of the Problem
A steady two-dimensional free convective flow in laminar form is considered over a vertical permeable circular cone of radius r with uniform heat flux. The cone is immersed in a viscous and incompressible fluid with temperature-dependent viscosity, for the numerical investigation. It is assumed that the surface heat flux of the cone is q w . Here, T ∞ is the ambient temperature of the fluid, and the considered geometry is illustrated in Figure 1.
the ambient temperature of the fluid, and the con 1. The equations governing the flow are [23][24][25]30,37] Energy equation : u ∂T ∂x The boundary conditions of Equations (1)-(3) are provided below [23]: where (u, v) are velocity components along the (x, y) axes, g is the acceleration due to gravity, ρ is the density, γ is the cone apex half-angle, k is the thermal conductivity, β is the coefficient of thermal expansion, µ(T) is the temperature-dependent viscosity of the fluid, where T denotes the temperature, and V is the transpiration velocity, which is positive for suction and negative for the injection of fluid through the cone. In the current research, only the suction case has been considered, and therefore V w has been considered to be positive throughout. Viscosity variations have multifarious forms. However, the one proposed by Charraudeau [38] is one of the most accepted ones and has been taken into account in this study. The expression could be presented as the following: where the suffix fil represents the film temperature of the fluid. If suction is non-existent in the model, the boundary layer acts as the free convection boundary layer near the leading edge, although improved downstream suction will be able to dominate the flow to a greater extent. As a result, the following transformations are required [23]: where ν ∞ (=µ ∞ /ρ) is the reference kinematic viscosity, Gr x is the local Grashof number, ξ is the dimensionless transpiration parameter, η is the pseudo-similarity variable, f (ξ,η) and θ(ξ,η) are the non-dimensional stream and temperature function, respectively, and ψ is the stream function defined by the following: Substituting (6) into Equations (1)-(5) and after some algebraic calculation, the following transformed equations are obtained: And the boundary conditions are presented hereby: where ε is the viscosity-variation parameter, and Pr is the Prandtl number defined as the following: The local non-similar partial differential equations presented in Equations (8) and (9) need to be solved after finding the boundary conditions (Equation (10)), where the latter could be obtained by the Keller box method, as verified in the literature [25]. The results are computed by considering uniform grids in the ξ-direction with 1001 grid-points considered in the 0 ≤ ξ ≤ 20. Further information on iterations and simulation strategies has been included in Section 3.1.
After calculating the respective values of f, θ, and their derivatives, the fundamental approach is to calculate the local skin-friction coefficient and the local Nu number from the following expressions: The results obtained by this method are presented in tabular form in Table 1 for different values of the viscosity variation parameter ε (=0.0, 0.5, 2.0) and Prandtl number Pr (=0.05, 0.1, 0.7). In the following sections, the solutions for small and large ξ have been presented, and are numerically valid and accurate in both the neighboring leading edge (small ξ) and downstream region (ξ). The value ξ for x is small near the leading edge or small V or both, and the series solution of Equations (8) and (9) can be found by applying the perturbation method considering x as a perturbation parameter. In order to establish such an objective, the functions f (ξ,η) and θ(ξ,η) in powers of ξ are expanded, and therefore: Equations (15)-(23) are solved pair-wise one after another. The solutions are obtained by combining the famous Runge-Kutta-Butcher [39] initial value solver with the Nachtseim-Swigert iteration scheme [40]. Thus, solutions are found for f i and θ i (i = 0, 1, 2) and their respective derivatives.
Once the values of f i and θ i for i = 0, 1, 2 and their derivatives are obtained, the local skinfriction coefficient and the local Nusselt number are calculated from the following expressions: The comprehensive values calculated from Equations (24) and (25) are presented in Table 1. The comparison also served as part of the validations of the present approach.

Solution for Large ξ
Attention has been given in this part to the behavior of the solution to the Equations (8) and (9) when ξ is considerably large. By the order of magnitude analysis of the various terms, ξθ f " in (8) and ξθ in (9) were found to be the largest. However, in both of their equations, both numerical terms need to be balanced mathematically. The balancing part is performed by assuming η to be a small parameter, which will eventually make η-derivatives larger. It is also important to determine the standard scaling approach considering θ = O(ξ − 1) as ξ→ ∞. After balancing the f , θ, and ξθ f " terms in (9), η = O(ξ − 1) and f = O(ξ − 4) as ξ→ ∞ would be found, which would also serve as a confirmation of the accuracy in the balancing. As a result, the following expressions are substituted: Substituting this transformation into Equations (8) and (9), we get the following equations: The corresponding boundary conditions are where primes describe differentiation with respect to η. Equations (27) and (28) are solved in terms of an inverse power of series in ξ.
The functions f (ξ, η) and θ(ξ, η) are expanded in the power series in the negative powers of ξ, considering that ξ is large, which yields the following: Now substituting the above expansions into Equations (27)- (29) and equating the coefficient of power of ξ 0 and ξ 1 , we get the following equations The local skin-friction and the local Nusselt number are as follows: The asymptotic solutions obtained from (31)- (36) in terms of local skin-friction and local Nusselt number and compared with the solution of the finite difference method in Table 1.

Computational Methods
The analyses have been conducted by tuning the non-dimensional empirical parameters at different stages with both fixed and variable Pr numbers. The computation was done by Fortran 90 [41]. The finite difference solutions were obtained by the Keller box scheme [31][32][33][34][35][36][37][38][39][40][41][42], which is considered to be one of the most accurate implicit finite difference techniques in computational mathematics. Refer to [37] for a detailed explanation of the Keller box method. In this method, the non-linear system of the partial differential equations governing the fluid flow is solved, and an assumption is made in terms of the functions and the derivatives to express the first-order equations. The derivatives in terms of the central differences are approximated in both co-ordinate directions. This part is particularly obtained by denoting the mesh points in the (ξ, η) plane by x i and η I , where i = 1, 2, 3, . . . , M and j = 1, 2, 3, . . . , N [37], followed by central difference approximations where the equations containing x are centred explicitly at (ξ i−1/2 , η j−1/2 ) and the rest at (ξ i , η j−1/2 ), where η j−1/2 = (η j + η j−1 )/2, for instance. It yields a batch of non-linear differential equations for the unknowns at xi in as a function of ξ i−1 . Then, those equations are subject to linearization by Newton's quasi-linearization method with a view to solving them by taking a block-tridiagonal algorithm into account. The initial iteration of the converged solution is considered to be at ξ = ξ i−1 . At the commencement, i.e., ξ = 0, the guess profiles for all the considered variables are provided, and then the Keller box method comes into the computation to solve the governing ordinary differential equations. By obtaining the lower stagnation point solution, the steps could be performed along with the boundary layer of the geometry [37]. The iteration is terminated once the target difference in velocity and temperature computation is reached, and in this study it was |δf i | ≤ 10 −6 , where i represents the number of iterations. It should be mentioned here that computations were not conducted by a uniform grid in the y-direction. However, a non-uniform grid was considered, which was then defined by η j = sinh(( j − 1)/p), with j = 1, 2, . . . , 301 and p = 100.

Overview of Numerical Analyses
At first, the validations were conducted by comparing the three different types of solutions considered in this study. A different type of 3D analyses was conducted to assess the Computation 2022, 10, 60 9 of 21 correlation coefficient, and good R 2 values were obtained overall. As part of the discussion, the streamlines and isothermal behavior have been presented. After that, the influence of parameters has been investigated on the viscosity, velocity, and temperature distributions. At the end of each analysis, a conclusion has been drawn outlining the significant and non-significant impact of specified parameters under different circumstances.

Validation and Comparison at Fixed Pr
To assess the accuracy of the approach, the typical finite difference solutions have been compared against the present approach with small and large suction parameters (ξ). To achieve such an objective, two different values of ε were assigned (0 and 2), and Pr = 0.7 was considered. In other words, the comparisons have been made in both the absence and presence of ε for the same type of fluid. The cross-validations were conducted by observing the values of two different parameters namely, the local skin-friction coefficient (C fx Gr x 1/5 ) and the local Nusselt number (Nu x Gr x −1/5 ). Table 1 contains the comprehensive comparisons of the aforementioned conditions. In general, there is hardly any big difference between the finite difference solutions and the asymptotic solutions. The slight deviation in the values could be attributed to the percentage error in the approach.
According to Table 1, as ξ increased, C fx Gr x 1/5 decreased, regardless of the ε values.
However, C fx Gr x 1/5 values were found to be lower in the presence of ε, thus confirming the role of the viscosity-variation parameter on the reduction of the local skin-friction coefficient. This behavior could be explained in terms of the viscosity and temperature of the fluid. As ε increased, the fluid temperature within the boundary layer increased, which led to increased viscosity. Therefore, C fx Gr x 1/5 decreased further. Similar patterns were observed for both finite difference and asymptotic solutions.
Meanwhile, Nu x Gr x −1/5 exhibited erratic trends. In general, Nu x Gr x −1/5 increased as ξ increased when ε = 0 was considered. However, as ε = 2 was assigned, Nu x Gr x −1/5 values had a different pattern. Initially, Nu x Gr x −1/5 kept on increasing as ξ increased when ε = 2, and the values of local Nu number were lower than those achieved at ε = 0. However, at ξ ≥ 2, the values of Nu x Gr x −1/5 exhibited values higher than Nu x Gr x −1/5 at ε = 0. Furthermore, surface analyses were conducted to understand the rational behavior as well as the accuracy of the approach. The purpose was to add one type of extra validation in the current approach to check the model predictions for a different type of fluid. Therefore, the local skin-friction coefficients and local Nu number were calculated at Pr = 0.7 and the numerical accuracy was assessed by observing the coefficient of correlation (R 2 ) at ε = 0 and ε = 2. The surface analyses were performed by the data analysis software Origin Pro, developed by OriginLab Corporation [42]. In the 3D analyses, solutions from the finite difference method and perturbation-asymptotic solutions were illustrated in the same frame to observe the analytical behavior as a function of the suction parameter (ξ). In general, an R 2 between 0.97 and 0.98 was obtained, which provides more confidence in the calculative approach. Figure 2 depicts the predicted surfaces for C fx Gr x 1/5 at different ξ in the absence of ε ( Figure 2a) and with an assigned value of ε ( Figure 2b). Most of the obtained points were found on the surface. In short, the marginal calculation error was noticed. On the other hand, a similar analysis was conducted in terms of Nu x Gr x −1/5 in Figure 3 under the same ε conditions. It is evident that Nu x Gr x −1/5 is increasing concurrently as ξ was increasing in Figure 3a,b, in line with the argument presented in Table 1, where Pr = 0.2. Furthermore, the calculated values from the asymptotic and perturbation solutions for small and large ξ were found to be closer than those from the finite difference approach. The missing values have not been included in the illustrations, and therefore the input sample size is not consistent in all figures. ure 3 under the same ε conditions. It is evident that NuxGrx −1/5 is increasing concurrently as ξ was increasing in Figure 3a,b, in line with the argument presented in Table 1, where Pr = 0.2. Furthermore, the calculated values from the asymptotic and perturbation solutions for small and large ξ were found to be closer than those from the finite difference approach. The missing values have not been included in the illustrations, and therefore the input sample size is not consistent in all figures.    increased concurrently. As ξ increased, the temperature of the fluid increased as mentioned earlier, and hence at ξ = 20 the values of C fx Gr x 1/5 were close to 0 due to the dominance of the suction parameter, whereas at the same suction parameter value, Nu x Gr x −1/5 was found to be the highest regardless of Pr values. However, it should be mentioned that as ξ kept increasing, the difference between the Nu x Gr x −1/5 of the fluid corresponding to Pr = 0.1 became wider than the former, with Pr = 0.05. It was observed that as Pr increased CfxGrx 1/5 values decreased, while NuxGrx −1/5 increased concurrently. As ξ increased, the temperature of the fluid increased as mentioned earlier, and hence at ξ = 20 the values of CfxGrx 1/5 were close to 0 due to the dominance of the suction parameter, whereas at the same suction parameter value, NuxGrx −1/5 was found to be the highest regardless of Pr values. However, it should be mentioned that as ξ kept increasing, the difference between the NuxGrx −1/5 of the fluid corresponding to Pr = 0.1 became wider than the former, with Pr = 0.05.

Comparison at Variable Pr
Since the temperature of the fluid kept increasing at increasing ξ values, the shear stress, i.e., the local skin-friction coefficient, decreased, and to balance the temperature difference, the rate of the heat transfer, i.e., the local Nu number, increased.  Since the temperature of the fluid kept increasing at increasing ξ values, the shear stress, i.e., the local skin-friction coefficient, decreased, and to balance the temperature difference, the rate of the heat transfer, i.e., the local Nu number, increased.

Development of Streamlines
This part of the discussion refers to fluids corresponding to two different types of Pr numbers (Pr = 0.05, 0.1) to depict the changes in the streamlines as a function of the suction parameter ξ. The streamlines illustrate the flow pattern within the momentum boundary layer. The changes were analyzed by observing the variations in the values of the flow rate (ψ). Two different viscosity-variation parameters (ε = 0.0, 0.5) were considered for each type of fluid to include the effect of ε on the streamlines. It was anticipated that the differences in the shapes of streamlines in both ε values would be quite marginal at low ξ(ξ < 4), due to the lower effect of the suction parameter. However, as ξ increases to a considerably larger value (ξ > 4), the differences between the shapes of streamlines of ε = 0.0 and ε = 0.5 increase and do not overlap. Furthermore, since ε is the viscosity-variation parameter, its absence will lead to the maximum flow rate, regardless of the Pr values. This behavior was expected due to the effect of viscosity on the fluid movement, which has been explained in the later sections, where the distribution of velocity profiles has been comprehensively discussed.
According to Figure 5, at lower ξ the streamlines corresponding to both ε values overlap, indicating the non-significant impact of the suction parameter on the fluid characteristics. A similar behavior was found for both Pr = 0.05 (Figure 5a) and Pr = 0.1 (Figure 5b). However, as ξ kept on increasing, the differences in the maximum values of the fluid flow rate of each segment corresponding to ε = 0.0 and ε = 0.5 were increasing gradually, which was mentioned in the first paragraph. This behavior could be also explained in terms of the effect of ε on the streamline's development. In the presence of the viscosity-variation parameter, the streamlines of the fluid seem to occupy less surface area due to the effect of ε, which eventually restricts the local maximum of the velocity distribution of the fluid. As per Figure 5a, while the maximum flow rate was found to be 42.90 in the absence of ε, the value plummeted to 13.11 (Figure 5b). On the other hand, the fluid corresponding to Pr = 0.1 exhibited a different trend in the pattern of streamlines. While the streamlines corresponding to the maximum flow rate were more aligned to the right side of the boundary layer, the maximum flow rate curve at Pr = 0.1 shifted to the center of the boundary, indicating the effect of Pr on the fluid characteristic. If Pr is increased, the momentum and thickness of the thermal boundary layer decreased, and hence the flow rate got reduced. of the fluid. As per Figure 5a, while the maximum flow rate was found to be 42.90 in the absence of ε, the value plummeted to 13.11 (Figure 5b). On the other hand, the fluid corresponding to Pr = 0.1 exhibited a different trend in the pattern of streamlines. While the streamlines corresponding to the maximum flow rate were more aligned to the right side of the boundary layer, the maximum flow rate curve at Pr = 0.1 shifted to the center of the boundary, indicating the effect of Pr on the fluid characteristic. If Pr is increased, the momentum and thickness of the thermal boundary layer decreased, and hence the flow rate got reduced.

Isothermal Behavior within the Boundary Layer
After investigating the behavior of the streamlines, a similar parametric effect was studied in terms of isotherms, illustrated in Figure 6. The isotherms depict the temperature distribution within the boundary layer. It was expected that the isothermal lines closer to the wall would be equal to the boundary temperature, which would lead to lines

Isothermal Behavior within the Boundary Layer
After investigating the behavior of the streamlines, a similar parametric effect was studied in terms of isotherms, illustrated in Figure 6. The isotherms depict the temperature distribution within the boundary layer. It was expected that the isothermal lines closer to the wall would be equal to the boundary temperature, which would lead to lines without any significant peak values. The isothermal line considerably far from the wall boundary would have peak values before plummeting gradually as the lines tended to come near the boundary layer.

R REVIEW
13 of 21 without any significant peak values. The isothermal line considerably far from the wall boundary would have peak values before plummeting gradually as the lines tended to come near the boundary layer. The isotherms appended in Figure 6 confirm the aforementioned statement. The isothermal lines far from the boundary layer could be seen with a peak before gradually starting to reduce. There is a slight deviation in the values between ε = 0.0 and ε = 0.5, where the presence of ε reasonably increased the isothermal lines due to the influence of the viscosity-variation. However, this was only visible when the lines were at an optimal distance from the thermal boundary layer, where the influence of the latter would be less The isotherms appended in Figure 6 confirm the aforementioned statement. The isothermal lines far from the boundary layer could be seen with a peak before gradually starting to reduce. There is a slight deviation in the values between ε = 0.0 and ε = 0.5, where the presence of ε reasonably increased the isothermal lines due to the influence of the viscosity-variation. However, this was only visible when the lines were at an optimal distance from the thermal boundary layer, where the influence of the latter would be less significant. A similar behavior was obtained in terms of streamlines as well where the presence of ε shifted the streamlines up compared to the streamlines corresponding to ε = 0.0. However, as the fluid motion decreases with an increase from Pr = 0.05 (Figure 6a) to Pr = 0.1 (Figure 6b), the peak values of isothermal lines decrease (visible ones) due to the reduction of the flow rate. The significant reduction in the maximum flow rate was also recorded in Figure 6b, and hence the consistency of the streamlines and isothermal analyses could be confirmed before heading towards the detailed parametric analyses in the following sections.
3.6. Impact of ε, η, and ξ Fluid Characteristics at Fixed Pr 3.6.1. Influence on Viscosity Distribution Figure 7 depicts the effects of both the viscosity-variation parameter (ε) and the suction parameter (ξ) on the dimensionless viscosity. The pseudo-similarity variable (η) was varied from 0 to a minimum of 6 (for ξ = 10) and a maximum of 10.0 (ξ = 1), depending on where the viscosity started to remain indifferent regardless of the η values. In Figure 7a,b, Pr = 0.1 was considered to maintain the consistency of the fluid characteristics.

FOR PEER REVIEW
14 of 21 = 10, the viscosity reached the unity at a lower η than that of ξ = 1, which was also anticipated.

Effect on Velocity Distribution (f')
In this part of the study, the changes in the non-dimensional velocity profiles were investigated by varying η and ε ∈ [0,2], at Pr = 0.1.
It was anticipated that, in the absence of ε, the velocity distribution would keep increasing as η increased and quickly reached the peak value, before η started to dominate, which would lead to a rapid decrease in the velocity (f'). As a consequence, at one specific value of η, the fluid would not exhibit any rapid mobility within the boundary. However, the local maximum of velocity would occur at a later stage of η if ε kept increasing. Furthermore, the static condition of fluid was also expected to be quicker for a higher ξ (in this case, ξ = 10), as a higher suction parameter demobilizes the fluid motion to a greater extent. Figure 8 illustrates the influence of the aforementioned parameters on the velocity distributions. As per Figure 8, as ε keeps increasing, the curves corresponding to the velocity exhibited peak values and later started to decrease significantly, heading to the static state as η increased. In general, it was found that as η increased, viscosity exhibited the local maximum values at the maximum ε, in this case ε = 2, as well as at η = 0. However, the local maximum viscosity was obtained at approximately 20.5 at ξ = 10, whereas it was found to be around 7.6 at ξ = 1. This behavior could be attributed to the fact that as the suction parameter increases, the viscosity exhibits the maximum values due to the absence of a pseudosimilarity variable which acts as a counterincentive to viscosity. On the other hand, ξ works as a catalyst to the viscosity at an increased suction parameter, and at η = 0, the temperature of the fluid within the respective boundary layer decreases, leading to the increased value of viscosity. Meanwhile, the curves corresponding to ε ∈ [0,2] in Figure 7a,b showed similar decreasing trends in viscosity. At ε = 0, there is no variation in the viscosity, and therefore η will always be 1, regardless of the values of other influential parameters. As ε began to increase, peak values started to increase from η = 0, before merging towards the unity, indicating an insignificant impact of η on the viscosity. Due to the dominance of ξ = 10, the viscosity reached the unity at a lower η than that of ξ = 1, which was also anticipated.

Effect on Velocity Distribution (f )
In this part of the study, the changes in the non-dimensional velocity profiles were investigated by varying η and ε ∈ [0,2], at Pr = 0.1.
It was anticipated that, in the absence of ε, the velocity distribution would keep increasing as η increased and quickly reached the peak value, before η started to dominate, which would lead to a rapid decrease in the velocity (f ). As a consequence, at one specific value of η, the fluid would not exhibit any rapid mobility within the boundary. However, the local maximum of velocity would occur at a later stage of η if ε kept increasing. Furthermore, the static condition of fluid was also expected to be quicker for a higher ξ (in this case, ξ = 10), as a higher suction parameter demobilizes the fluid motion to a greater extent. Figure 8 illustrates the influence of the aforementioned parameters on the velocity distributions. As per Figure 8, as ε keeps increasing, the curves corresponding to the velocity exhibited peak values and later started to decrease significantly, heading to the static state as η increased.

Effect on Velocity Distribution (f')
In this part of the study, the changes in the non-dimensional velocity profiles were investigated by varying η and ε ∈ [0,2], at Pr = 0.1.
It was anticipated that, in the absence of ε, the velocity distribution would keep increasing as η increased and quickly reached the peak value, before η started to dominate, which would lead to a rapid decrease in the velocity (f'). As a consequence, at one specific value of η, the fluid would not exhibit any rapid mobility within the boundary. However, the local maximum of velocity would occur at a later stage of η if ε kept increasing. Furthermore, the static condition of fluid was also expected to be quicker for a higher ξ (in this case, ξ = 10), as a higher suction parameter demobilizes the fluid motion to a greater extent. Figure 8 illustrates the influence of the aforementioned parameters on the velocity distributions. As per Figure 8, as ε keeps increasing, the curves corresponding to the velocity exhibited peak values and later started to decrease significantly, heading to the static state as η increased.  At ε = 0 and ξ = 1 (Figure 8a), the velocity reached its maximum value of 0.955 at η = 1, whereas it was 0.85 at ε = 0.5. The lowest peak was recorded at the maximum ε values considered in this section (ε = 2), and due to the dominance of η and ε, the fluid had the lowest peak. In any case, after reaching the peak, velocity started to decrease significantly before reaching the completely static state at η = 10. On the other hand, at ξ = 10 (Figure 8b), similar trends were observed for the velocity at different ε and η, but due to the 10-time-increasing high suction parametric value, the peak velocity values were significantly lower, which indicates the superiority of the high ξ values. It is also expected that further increases in ξ will lead to further reduced peak values, and, at one stage, the fluid will remain static.

Changes in Temperature Distribution
Studying the changes in the temperature was crucial to the current research, since it provided more information on the effect of ε and η. To maintain consistency in the presentation, the changes in the non-dimensional temperature were studied by varying η and ε ∈ [0,2], at Pr = 0.1. Figure 9 reports the consistent decreasing trend of temperature as η increased. However, the peak θ was recorded at the highest ε (in this case ε = 2) and lowest pseudo-similarity variable (η = 0), where ξ was unity. This behavior could be well-explained in the light of a conceptual understanding of the effect of the viscosity-variation parameter on the temperature. In the absence of η and ε, the temperature exhibited the lowest peak value of 2.9, whereas at ε = 2 the value of dimensionless temperature was found to be 3.35. The differences in the θ outline the impact of the viscosity-variation parameter on the temperature. However, as η started to have non-zero input, the temperature kept on decreasing rapidly and it was close to 0 as η reached maximum input values corresponding to the case studies (η = 10 for Figure 9a, η = 6 for Figure 9b). The lowest temperature values could also be linked with the velocity profiles obtained in Figure 8, which showed that in similar η values, the velocity reached 0, thus almost demobilizing the fluid within the boundary layer. Since there was no velocity or marginal velocity recorded at η = 10 in Figure 8a and η = 6 in Figure 8b, this backs up the concept that temperature will be close to 0 or equal to 0 as the fluid remains static. In addition, the influential impact of the suction parameter ξ largely reduced the temperature of the fluid, and the impact of ε was quite marginal, indicating the dominance of ξ (Figure 9b). While the curves seemed to overlap in Figure 9b, the magnified frame demonstrated the marginal differences at the narrowed axes. ever, the peak θ was recorded at the highest ε (in this case ε = 2) and lowest pseudosimilarity variable (η = 0), where ξ was unity. This behavior could be well-explained in the light of a conceptual understanding of the effect of the viscosity-variation parameter on the temperature. In the absence of η and ε, the temperature exhibited the lowest peak value of 2.9, whereas at ε = 2 the value of dimensionless temperature was found to be 3.35. The differences in the θ outline the impact of the viscosity-variation parameter on the temperature. However, as η started to have non-zero input, the temperature kept on decreasing rapidly and it was close to 0 as η reached maximum input values corresponding to the case studies (η = 10 for Figure 9a, η = 6 for Figure 9b). The lowest temperature values could also be linked with the velocity profiles obtained in Figure 8, which showed that in similar η values, the velocity reached 0, thus almost demobilizing the fluid within the boundary layer. Since there was no velocity or marginal velocity recorded at η = 10 in Figure 8a and η = 6 in Figure 8b, this backs up the concept that temperature will be close to 0 or equal to 0 as the fluid remains static. In addition, the influential impact of the suction parameter ξ largely reduced the temperature of the fluid, and the impact of ε was quite marginal, indicating the dominance of ξ (Figure 9b). While the curves seemed to overlap in Figure 9b, the magnified frame demonstrated the marginal differences at the narrowed axes.  In general, the Pr number distinguishes different types of fluids. For example, Pr = 0.01 refers to sodium, Pr between 0.02 to 0.03 refers to mercury. The purpose of different Pr tests in this study is to explore the uniformity of the numerical solution. Figure 10 shows the changes in the non-dimensional viscosity values as a function of Pr. The value of ε was kept constant at 1, whereas two different suction parameters were appended in the numerical simulations. In addition, η ∈ [0,10] was also considered. It was observed that the fluid corresponding to Pr = 0.01 exhibited the highest peak value when η = 0, indicating the absence of η or the fact that the low viscosity-variation parameter increased the viscosity value. As Pr got increased, the peak values got reduced as well. While the local maximum viscosity of the fluid of Pr = 0.01 was found to be 9.45 at η = 0, the fluid corresponding to Pr = 0.1 exhibited a local maximum viscosity of 4.2. However, as η kept on increasing, the viscosity kept decreasing due to the dominance of η on the viscosity values. However, these phenomena were obtained at a low suction parametric value ξ = 1, as presented in Figure 10a. As ξ increased to 5, the peak values of each fluid at different Pr exhibited a rapid increase in the local maxima. The fluid corresponding to Pr = 0.01 demonstrated a peak viscosity value of approximately 40, which was 9.45 when ξ = 1 was considered. This 323.28% increase in the value could be attributed to the increase in ξ, which eventually increased the viscosity value in the absence of η. The rest of the fluids with further increased Pr also demonstrated similar characteristics, as depicted in Figure 10b. as η kept on increasing, the viscosity kept decreasing due to the dominance of η on the viscosity values. However, these phenomena were obtained at a low suction parametric value ξ = 1, as presented in Figure 10a. As ξ increased to 5, the peak values of each fluid at different Pr exhibited a rapid increase in the local maxima. The fluid corresponding to Pr = 0.01 demonstrated a peak viscosity value of approximately 40, which was 9.45 when ξ = 1 was considered. This 323.28% increase in the value could be attributed to the increase in ξ, which eventually increased the viscosity value in the absence of η. The rest of the fluids with further increased Pr also demonstrated similar characteristics, as depicted in Figure 10b.

Changes in the Velocity Distribution
Velocity distribution corresponding to different Pr implies the influence of Pr on different types of fluids, as other parameters, such as η and ξ, are varied. The viscosity-variation parameter ε was fixed at 1. In general, the patterns of velocity distribution had similarities with those obtained in Figure 11, which also indicates consistency in the current approach. As η started to increase from 0, the velocity started to increase rapidly before hitting the peak, followed by a sharp decrease in a fluid motion. However, within the η ∈ [0,10] boundary, only fluids corresponding to Pr = 0.01, 0.02 remained mobile, although the values of f' were significantly lower at η = 10, as shown in Figure 11a. On the other hand, the rest of the fluids headed towards the static condition within the boundary layer at η = 10, thus establishing the preeminence of η over the fluid mobility. A similar behavior was observed at an increased ξ = 5 (Figure 8b), but the influence on f' was not robust. As ξ increased to 5 from 1 in Figure 11b, the local maxima decreased. While at ξ = 1, Pr = 0.02,

Changes in the Velocity Distribution
Velocity distribution corresponding to different Pr implies the influence of Pr on different types of fluids, as other parameters, such as η and ξ, are varied. The viscosityvariation parameter ε was fixed at 1. In general, the patterns of velocity distribution had similarities with those obtained in Figure 11, which also indicates consistency in the current approach. As η started to increase from 0, the velocity started to increase rapidly before hitting the peak, followed by a sharp decrease in a fluid motion. However, within the η ∈ [0,10] boundary, only fluids corresponding to Pr = 0.01, 0.02 remained mobile, although the values of f were significantly lower at η = 10, as shown in Figure 11a. On the other hand, the rest of the fluids headed towards the static condition within the boundary layer at η = 10, thus establishing the preeminence of η over the fluid mobility. A similar behavior was observed at an increased ξ = 5 (Figure 8b), but the influence on f was not robust. As ξ increased to 5 from 1 in Figure 11b, the local maxima decreased. While at ξ = 1, Pr = 0.02, the local maximum was recorded to be approximately 1.4, and the value was found to be around 0.85 when ξ = 5 was considered. It should be mentioned here that while ξ attenuated the peak values of the velocity, it also influenced the local maxima of fluids with lower Pr values to occur at higher η. For example, fluid corresponding to Pr = 0.02 had a peak at η ≈ 2 at ξ = 1 (Figure 11a), whereas it was found to be at η ≈ 3 when ξ = 5 (Figure 11b).

Effect on Temperature Distribution
In the final part of sensitivity analyses, the effect of dimensionless parameters η and ξ has been investigated on the temperature distribution as Pr varied from 0.01 to 0.1 under no viscosity variation, i.e., ε was constant 1 throughout this analysis. In general, it was expected that η would have more influence on the temperature distribution with varied Pr numbers, as fluids' temperature would likely to be more influenced by Pr values, as it would characterize the type of fluid, and η would be able to dominate, since the viscosityvariation parameter was kept constant. The suction parameter ξ was expected to have a marginal influence on the temperature distribution of fluids with low Pr (in this case 0.01 and 0.02) and a pronounced effect on the rest. In any case, the influence of ξ will not be significantly large on temperature distribution. the local maximum was recorded to be approximately 1.4, and the value was found to be around 0.85 when ξ = 5 was considered. It should be mentioned here that while ξ attenuated the peak values of the velocity, it also influenced the local maxima of fluids with lower Pr values to occur at higher η. For example, fluid corresponding to Pr = 0.02 had a peak at η ≈ 2 at ξ = 1 (Figure 11a), whereas it was found to be at η ≈ 3 when ξ = 5 ( Figure  11b).

Effect on Temperature Distribution
In the final part of sensitivity analyses, the effect of dimensionless parameters η and ξ has been investigated on the temperature distribution as Pr varied from 0.01 to 0.1 under no viscosity variation, i.e., ε was constant 1 throughout this analysis. In general, it was expected that η would have more influence on the temperature distribution with varied Pr numbers, as fluids' temperature would likely to be more influenced by Pr values, as it would characterize the type of fluid, and η would be able to dominate, since the viscosityvariation parameter was kept constant. The suction parameter ξ was expected to have a marginal influence on the temperature distribution of fluids with low Pr (in this case 0.01 and 0.02) and a pronounced effect on the rest. In any case, the influence of ξ will not be significantly large on temperature distribution. Figure 12 supports the aforementioned claims. For instance, the fluid corresponding to Pr = 0.01 had the highest peak temperature value, regardless of the ξ values. However, the peak was only visible at η = 0 in both cases (ξ = 1 or ξ = 5). As η increased gradually, the temperature value θ kept decreasing, indicating a reduction in fluid temperature within the boundary layer. However, fluids with Pr ≥ 0.05 exhibited a θ close to 0 as η = 10 was considered in the input. This behavior was also observed in both ξ values. At Pr = 0.01, θ marginally decreased from 7.9 (Figure 12a) to 7.6 ( Figure 12b) as ξ increased from 1 to 5 when η was non-existent in the computation. A similar attribute was noticed at Pr = 0.02 as well. However, at Pr = 0.1, the plummeting rate was significantly larger. While θ = 3.6 was observed for P r= 0.1 under no η condition (Figure 12a), the value was found to be decreased to 1.8 when ξ = 5 was considered. This decrease could be attributed to the thermal properties of the fluid. Increased Pr values imply the fluids which have low or decreasing thermal conductivity. Therefore, as the Pr increased, the fluid was not able to conduct heat significantly, which eventually decreases the thickness of the boundary layer.  Figure 12 supports the aforementioned claims. For instance, the fluid corresponding to Pr = 0.01 had the highest peak temperature value, regardless of the ξ values. However, the peak was only visible at η = 0 in both cases (ξ = 1 or ξ = 5). As η increased gradually, the temperature value θ kept decreasing, indicating a reduction in fluid temperature within the boundary layer. However, fluids with Pr ≥ 0.05 exhibited a θ close to 0 as η = 10 was considered in the input. This behavior was also observed in both ξ values. At Pr = 0.01, θ marginally decreased from 7.9 (Figure 12a) to 7.6 ( Figure 12b) as ξ increased from 1 to 5 when η was non-existent in the computation. A similar attribute was noticed at Pr = 0.02 as well. However, at Pr = 0.1, the plummeting rate was significantly larger. While θ = 3.6 was observed for P r= 0.1 under no η condition (Figure 12a), the value was found to be decreased to 1.8 when ξ = 5 was considered. This decrease could be attributed to the thermal properties of the fluid. Increased Pr values imply the fluids which have low or decreasing thermal conductivity. Therefore, as the Pr increased, the fluid was not able to conduct heat significantly, which eventually decreases the thickness of the boundary layer.

Conclusions
The present study discusses the effect of viscosity-variation in the natural convective

Conclusions
The present study discusses the effect of viscosity-variation in the natural convective flow inside a permeable vertical cone with a uniform surface heat flux. The numerical analyses performed in this study will provide elaborate explanations on the novelty of considering temperature-dependent viscosity to aid in different environmental and agricultural applications targeting uniform heat flux of a porous medium such as soil. Soil contains different types of fluids and may be represented in different geometries to explore all possible scenarios that occurred by either natural or human-made phenomena. Different types of fluids have been assigned by varying Pr numbers. Later, further investigations have been conducted in different parametric studies, which could be subcategorized into two major categories: (i) fixed Prandtl number (Pr) and (ii) variable Prandtl number (Pr). The numerical model was developed by transforming the governing boundary layer equations into a dimensionless form. In addition, the local non-similarity equations, obtained as a result of a non-linear system of partial differential equations' reduction, were solved by three different solution techniques. The validations were conducted at both fixed Pr and variable Pr to assess the accuracy of the model. The model assessment was conducted by comparing the finite difference solution with the asymptotic solutions for large and small ξ values, and good agreements were established. The summary of the present study is mentioned in the following:

•
Increasing the suction parameter (ξ) leads to decreasing shear stress (local skin-friction coefficient) and increasing the rate of heat transfer (local Nusselt number). The increasing and decreasing characteristics could be attributed to the temperature difference of the fluid within the boundary layer, which requires balancing the physical difference. • As the viscosity-variation parameter (ε) increases, the local skin-friction coefficient decreases concurrently due to the effect of viscosity. On the other hand, increasing the rate of heat transfer as a function of ε does not remain consistent due to the effect of ξ. A small ξ rate of heat transfer decreases as ε increases, but the opposite behavior is observed with a large ξ, which indicates the superiority of the suction parameter over viscosity.

•
If Pr increases, the local skin-friction coefficient decreases, and the rate of heat transfer increases rapidly due to the changes in the physical characteristics of the fluid and its viscosity.

•
At ε = 0, as η increases, viscosity decreases rapidly and heads towards 0 due to the dominance of η. However, at fixed ε and variable Pr, the curves corresponding to viscosity values exhibit the lowest value at a much later stage, indicating the high viscous characteristics of fluids with low Pr number. • At any value of ε, at η ∈ [0,1], velocity increases the maximum at one point and then sharply plummets towards static condition at η > 6. However, a large suction parameter (ξ = 10) significantly lowers the peak values regardless of ε or types of fluids. • An increased value of ε leads to the highest local maximum of temperature distribution in the absence of η and at a fixed Pr number. However, at η = 0, the temperature starts to decrease gradually. The large suction parameter (ξ = 10) also suggestively lowers the peak values regardless of ε or types of fluids. Meanwhile, at the variable Pr number, the local maxima of temperature distribution get marginally affected at a low Pr number (<0.05). Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Data availability on request.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Nomenclature
C

Specific heat C f
Skin-friction