Modeling and Analysis of 3D Magnetized Radiative Williamson Nanofluid Flow Past a Riga Plate with Forchheimer-Darcy Porous Medium and Cattaneo-Christov Double Diffusion

Non-Newtonian fluids are liquids that do not exhibit a direct linear relationship between shear stress and shear rate, unlike Newtonian fluids. This distinctive behaviour makes them fundamentally different from Newtonian fluids such as air and water, whose viscosity remains constant under given conditions. Non-Newtonian fluids display a wide range of complex flow characteristics and are classified based on their response to shear stress or shear rate, including shear-thinning, Bingham plastic, shear-thickening, thixotropic, viscoelastic, and rheopectic fluids. Determining their rheological properties is essential in several industries and applications, such as drilling operations, materials research, pharmaceuticals, cosmetics, and food processing. Traditional Navier-Stokes relations often fail to capture the complex behaviour of these fluids. Consequently, various mathematical models have been developed to describe their thermophysical and rheological characteristics, including the Carreau model [1], Maxwell fluid model [2], Cross model [3], Tiwari-Das model [4], Casson fluid model [5], Williamson model [6], among others. Sarda et al. [7] investigated heat and mass transfer in a non-Newtonian nanofluid through porous media considering thermophoresis and Brownian motion. Kumar et al. [8] carried out a computational study of 2D fluid flow with heat and mass transfer in the natural convection of a micropolar fluid subjected to magnetohydrodynamic (MHD) effects over a deformable, porous, stretched surface. Khan et al. [9] examined the influence of slip boundary conditions, modified Fourier and Fick laws, and Joule heating on the flow of a second-order fluid. Numerous researchers have explored the behavior of non-Newtonian fluids over stretching surfaces, as discussed in references [10], [11], [12], [13], [14], [15].

Nanofluids represent an advanced approach to heat transfer, consisting of a conventional base fluid into which nanoparticles are dispersed. They have been widely used in enzyme biosensors, microreactors, bio-separation systems, micro heat pipes, microchannel heat sinks, and various other devices. Over the past two decades, two-component nanofluid flows have attracted significant research attention because of their enhanced heat transfer capabilities. In many engineering and industrial applications, effective thermal management—reincluding cooling, heating, thermal efficiency, and energy savings-is of critical importance. Common base fluids such as hydrocarbons, water, ethylene glycol, oil, and mercury contribute substantially to these thermal processes. The addition of nanoparticles (1–100 nm) to a conventional fluid increases its thermal conductivity, forming what are known as nanofluids (Choi and Eastman [16]). Building on these developments, Buongiorno [17] proposed a comprehensive mathematical model that incorporates thermophoresis and Brownian motion effects. Water, due to its accessibility and favourable properties, is frequently used as a solvent in industrial operations. Gireesha et al. [18] examined the influence of thermal stratification on boundary-layer flow over an expanding sheet. Salleh et al. [19], using the Buongiorno nanofluid framework, theoretically investigated the roles of heat generation/absorption and chemical reactions in laminar boundary-layer nanofluid flow. Waleig et al. [20] analysed the influence of heat generation and chemical reactions on unsteady magnetohydrodynamic (MHD) radiative flow of a Williamson nanofluid over a permeable stretching surface. Waqas et al. [21] studied the impact of temperature and MHD effects on nanofluid flow in parallel-plate rotating systems. Sumathi et al. [22] explored the effects of key parameters on MHD nanofluid flow over an inclined stretching sheet. Numerous other researchers have also contributed significantly to nanofluid studies [23], [24], [25], [26], [27], [28], [29], [30], [31], [32].

Most significant fluid flows in chemical, manufacturing, biological, and industrial processes occur in the presence of porous or permeable media. Examples include petroleum purification, water penetration in storage tanks, material filtration, and even blood circulation. To study such systems effectively, it is important to investigate Darcy-Forchheimer fluid flow. While Darcy’s law provides a foundational description of flow through porous media, its applicability becomes limited at higher velocities. In many practical situations, especially when the porous region has low porosity or when the flow rate near the wall is high, Darcy’s linear relationship fails to describe the flow accurately. The Darcy-Forchheimer model incorporates both viscous and inertial effects, offering a more realistic representation of fluid velocity and heat transfer through porous structures. Considering exponential stretching, Ahmad et al. [33] examined a 3D couple-stress MHD nanofluid flow using the Darcy-Forchheimer model and reported that fluid acceleration increases with a rise in the inertia parameter. Awais et al. [34] investigated the influence of a Darcy-Forchheimer porous medium on the flow, temperature, and mass transfer characteristics of a 2D Eyring-Powell fluid. Influence of inertia “Forchheimer term” and macroscopic shear “Brinkman term” are pointedly considered by Elkady et al. [35]. Zolotukhin et al. [36] focused on focus of this article is on accounting for non-Darcy flow when inertial forces become more significant and the relationship between pressure gradient and seepage velocity becomes nonlinear.

A 3D flow model is a vital component of fluid flow analysis because it captures the full flow behavior observed in real-life applications, where variations occur in all three spatial directions. Unlike 1D or 2D analyses, a 3D analysis can represent detailed spatial gradients in velocity, pressure, and temperature, as well as important features such as secondary flows and certain turbulence effects that lower-dimensional models cannot capture. This level of detail is essential for making reliable predictions and guiding design improvements in fields such as aerodynamics (e.g., aircraft wing design), hydrodynamics (e.g., ship hull performance), and various engineering systems including turbines, motors, and heat exchangers. A 3D model provides more accurate estimates of forces, drag, lift, heat transfer, and energy losses, ultimately contributing to safer, more efficient, and cost-effective designs across aerospace, mechanical, and energy applications. Abdul Hakeem et al. [37] analyzed the effect of nonlinear radiation on nanofluid flow, Akolade and Tijani [38] compared two non-Newtonian nanofluids flowing over a Riga plate, and Jusoh et al. [39] reported the characteristics of 3D nanofluid flow over a permeable sheet.

An electromagnetic actuator that consists of the aggregation of adjusted magnets along with the spanwise arranged layout of electrodes with alternation mounting on a flat surface is indicated as a Riga plate. This actuator is responsible for the generation of exponentially decreasing Lorentz force within the boundary layer. Khatun et al. [40] performed on the numerical investigation of electro-magnetohydrodynamic (EMHD) radiating fluid flow nature along an infinitely long vertical Riga plate with suction in a rotating system. Reddy et al. [41] and Zaydan et al. [42] and focused on this subject.

The precise manipulation of heat conduction remains an intricate and pressing challenge. Historically, the classical Fourier’s law, grounded in experimental evidence, has been the prevailing model delineating the heat conduction process. However, a crucial limitation of Fourier’s law is its assumption of infinite heat conduction velocity, rendering it unsuitable for numerous applications. To bridge this gap, scholars have presented various extended models to portray heat conduction, notably Abro et al. [43], Benenti et al. [44], and Chen [45].

To the best of our knowledge, previous studies examined MHD, porous media, or Cattaneo-Christov flux separately (or in pairs) for various non-Newtonian models. However, the simultaneous coupling of (a) Williamson type shear-thinning rheology, (b) Cattaneo-Christov thermal and mass relaxation (non-Fourier effects), (c) Darcy-Forchheimer porous inertia, and (d) Riga surface electromagnetic forcing leads to non-trivial interactions: relaxation terms modify near-wall temperature and nanoparticle boundary layers (thus changing the effective Nusselt and Sherwood numbers), while porous inertial terms alter momentum redistribution that interacts with Lorentz forcing from the Riga plate. These cross-effects have not been catalogued collectively for a 3D flow; doing so elucidates parameter regimes where MHD actuation can be used to compensate for porous or rheological losses-a practical design insight for microfluidic cooling and magnetically assisted drug delivery. Table 1 shows the innovation of the existing study.

The coupled model developed in this study offers clear engineering value and practical applicability. It provides design guidance for MHD-assisted cooling and micro-scale heat exchangers by showing how magnetic/Riga actuation can enhance surface heat removal even in porous or shear-thinning fluids. It also helps optimize microfluidic and lab-on-a-chip systems, where the interaction between Darcy-Forchheimer drag and Lorentz forcing identifies conditions that restore flow rates and maintain required thermal gradients. In biomedical applications, the combined effects of Cattaneo-Christov relaxation and Buongiorno nanoparticle transport improve predictions of nanoparticle distribution and heating, which is important for magnetic hyperthermia and drug-delivery processes. Overall, the novelty of this work lies in integrating multiple advanced physical mechanisms within a 3D framework and translating their combined influence into practical, quantitative design guidelines for engineers working with MHD cooling devices, porous microreactors, and magnetically controlled biomedical flows.

The key objectives of this study are as follows (Figure 1).

i. To develop a 3D model for Williamson nanoliquid flow over a Riga surface with the Cattaneo-Christov structure by incorporating the Buongiorno model, and

ii. Transform the governing PDEs into a system of ODEs using similarity variables.

iii. To solve the system using Mathematica’s NDSolve technique, and

iv. Analyze the influence of various physical parameters on flow, heat, and mass transfer using graphical and tabular illustrations.

The results, presented through 2D and 3D plots, illustrate the variations in skin friction, and heat and mass transfer rates. A comparative analysis was performed to validate the accuracy of the numerical solution. This study offers operational insights that are valuable for engineering and industrial applications involving non-Newtonian nanofluids in porous materials. The findings are also relevant to biomedical engineering-for example, in targeted drug delivery and magnetic nanoparticle hyperthermia-where precise control of flow and heat distribution is critical. Additionally, chemical processing, microfluidic device design, and filtration technologies benefit from the use of magnetically active porous substrates.

The main goal of this acquisition is to scrutinize the incompressible, steady 3D radiative Williamson nanofluid flow over a Riga surface that is placed at $z=0$ and is stretching with the velocities $u=u_w(x)=a x$ and $v=v_w(x)=b x$ in $x$ and $y$ orientations respectively. Let $T_w$ and $C_w$ denote the surface temperature and concentration, respectively, which are invariably greater than the ambient temperature and concentration, symbolized by $T_{\infty}$ and $C_{\infty}$. The heat transfer and concentration features of the material were examined using the Cattaneo-Christov model. The effects of heat absorption/generation, thermal radiation, and permeable medium were also considered. The formulation properties of the nanofluids were the thermophoresis and Brownian motion.A convective heating process is responsible for controlling the surface temperature and, is characterized by the heat transport coefficient $h_f$. A physical model of the flow problem is shown in Figure ~\ref{fig2}.

Nadeem and Hussain [50] provided the mathematical structure of the Cauchy stress tensor for the Williamson fluid as follows:

$\boldsymbol{C}_s=-p \boldsymbol{I}+\bar{\tau}$, $\bar{\tau}=\left[\mu_{\infty}+\frac{\left(\mu_0-\mu_{\infty}\right)}{1-\Gamma \gamma^*}\right] A_1$, where $A_1=(\operatorname{grad} \boldsymbol{V})+(\operatorname{grad} \boldsymbol{V})^{T r}$.

From Nadeem and Hussian [50] we have: $\bar{\tau}=\mu_0\left[ 1+\Gamma \gamma^*\right] A_1$.

The governing equations of the flow are [39, 40]:

$ \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial \omega}{\partial z}=0 $

$ \left[u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}+\omega \frac{\partial u}{\partial z}\right]=v \frac{\partial^2 u}{\partial z^2}+\sqrt{2} v \Gamma \frac{\partial u}{\partial z} \cdot \frac{\partial^2 u}{\partial z^2}-\frac{v}{k_1} u-\frac{C_b}{x \sqrt{k_1}} u^2+\frac{\pi J_0 M_0 \exp \left[-\frac{\pi}{r^*} z\right]}{8 \rho} $

$ \left[u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+\omega \frac{\partial v}{\partial z}\right]=v \frac{\partial^2 v}{\partial z^2}+\sqrt{2} v \Gamma \frac{\partial v}{\partial z} \cdot \frac{\partial^2 v}{\partial z^2}-\frac{v}{k_1} v-\frac{C_b}{x \sqrt{k_1}} v^2 $

$ \begin{aligned} \left(u \frac{\partial T}{\partial x}+v \frac{\partial T}{\partial y}+\omega \frac{\partial T}{\partial z}\right)=\alpha_m & \frac{\partial^2 T}{\partial z^2}+\tau\left[D_B \frac{\partial C}{\partial z} \frac{\partial T}{\partial z}+\frac{D_T}{T_{\infty}}\left(\frac{\partial T}{\partial z}\right)^2\right] \\ & -\lambda_E\left[\frac{\partial^2 T}{\partial x^2}\left(u^2\right)+\frac{\partial^2 T}{\partial y^2}\left(v^2\right)+\frac{\partial^2 T}{\partial z^2}\left(\omega^2\right)+\frac{\partial^2 T}{\partial x \partial y} 2 u v\right. \\ & +\frac{\partial^2 T}{\partial y \partial z}(2 v \omega)+\frac{\partial^2 T}{\partial x \partial z}(2 u \omega)+\left(\frac{\partial u}{\partial x} u+\frac{\partial u}{\partial y} v+\frac{\partial u}{\partial z} \omega\right) \frac{\partial T}{\partial x} \\ & \left.+\left(u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+\omega \frac{\partial v}{\partial z}\right) \frac{\partial T}{\partial y}+\left(u \frac{\partial \omega}{\partial x}+v \frac{\partial \omega}{\partial y}+\omega \frac{\partial \omega}{\partial z}\right) \frac{\partial T}{\partial z}\right] \\ & +\frac{Q_0}{\rho c_p}\left(T-T_{\infty}\right)-\frac{1}{\left(\rho c_p\right)} \frac{\partial q_r}{\partial z} \end{aligned} $

$ \begin{aligned} u \frac{\partial C}{\partial x}+v \frac{\partial C}{\partial y}+\omega \frac{\partial C}{\partial z}= & D_B \frac{\partial^2 C}{\partial z^2}+\frac{D_T}{T_{\infty}} \frac{\partial^2 T}{\partial z^2}-\lambda_C\left[\frac{\partial^2 C}{\partial x^2}\left(u^2\right)+\frac{\partial^2 C}{\partial y^2}\left(v^2\right)+\frac{\partial^2 C}{\partial z^2}\left(\omega^2\right)\right. \\ & +\frac{\partial^2 C}{\partial x \partial y}(2 u v)+\frac{\partial^2 C}{\partial y \partial z}(2 v \omega)+\frac{\partial^2 C}{\partial x \partial z}(2 u \omega)+\left(\frac{\partial u}{\partial x} u+\frac{\partial u}{\partial y} v+\frac{\partial u}{\partial z} \omega\right) \frac{\partial C}{\partial x} \\ & \left.+\left(u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+\omega \frac{\partial v}{\partial z}\right) \frac{\partial C}{\partial y}+\left(u \frac{\partial \omega}{\partial x}+v \frac{\partial \omega}{\partial y}+\omega \frac{\partial \omega}{\partial z}\right) \frac{\partial C}{\partial z}\right] \end{aligned} $

The boundary constraints are:

$ \left.\begin{array}{c} u=a x, v=b y, \omega=0,-k \frac{\partial T}{\partial z}=h_f\left[T_w-T\right], C=C_w, \text { at } z=0 \\ u \rightarrow 0, \quad v \rightarrow 0, \quad T \rightarrow T_{\infty}, \quad C \rightarrow C_{\infty}, \quad \text { as } z \rightarrow \infty \end{array}\right\} $

Following Rosseland approximation, the radiative heat flux is: $q_r=-\frac{4 \sigma^* \partial T^4}{3 k^* \partial z}$.

Furthermore, $T^4$ is simplified to a linear function of temperature because the difference between the ambient temperature and fluid temperature within the flow is very small. By expanding $T^4$ in Taylor series about $T_{\infty}$, if we ignore higher-order terms, we obtain $T^4 \cong 4 T_{\infty}^3 T-3 T_{\infty}^4$.

Now the similarity transformations are:

$ \left.\begin{array}{l} u=a x f^{\prime}(\zeta), \zeta=\sqrt{\left(\frac{a}{v}\right)} z, v=a y g^{\prime}(\zeta), \omega=-\sqrt{a v}(f(\zeta)+g(\zeta)), \theta(\zeta)=\frac{T-T_{\infty}}{T_w-T_{\infty}}, \\ \phi(\zeta)=\frac{C-C_{\infty}}{C_w-C_{\infty}}. \end{array}\right\} $

Using Eq. (7), Eq. (2), and (5) can be interpreted as:

$ f^{\prime \prime \prime}+W f^{\prime \prime} f^{\prime \prime \prime}-f^{\prime 2}+[f+g] f^{\prime \prime}-\lambda f^{\prime}-F r f^{\prime 2}+H a e^{-\delta \zeta}=0 $

$ g^{\prime \prime \prime}+W g^{\prime \prime} g^{\prime \prime \prime}-g^{\prime 2}+[f+g] g^{\prime \prime}-\lambda g^{\prime}-F r g^{\prime 2}=0 $

$ \begin{aligned} & \frac{1}{P r}\left(1+\frac{4 R}{3}\right) \theta^{\prime \prime}+[f+g] \theta^{\prime}-\Gamma_T\left[[f+g]^2 \theta^{\prime \prime}+[f+g]\left[f^{\prime}+g^{\prime}\right] \theta^{\prime}\right]+N b \theta^{\prime} \phi^{\prime}+\\ & N t \theta^{\prime 2}+Q \theta=0 \end{aligned} $

$ \phi^{\prime \prime}+L e[f+g] \phi^{\prime}-L e \Gamma_c\left[[f+g]^2 \phi^{\prime \prime}+[f+g]\left[f^{\prime}+g^{\prime}\right] \phi^{\prime}\right]+\frac{N t}{N b} \theta^{\prime \prime}=0 $

Now the boundary constraints are

$ \left.\begin{array}{l} f(0)=0, g(0)=0, \phi(0)=1, f^{\prime}(0)=1, g^{\prime}(0)=\varepsilon, \theta^{\prime}(0)=-B i(1-\theta(0))\\ f(\infty)=g(\infty)=\theta(\infty)=\phi(\infty)=0 . \end{array}\right\} $

where,

$\begin{aligned} & \lambda=\frac{v}{k_1 a}, F r=\frac{C_b}{\sqrt{k_1}}, H a=\frac{\pi j_0 M_0}{8 a^3 x \rho}, \delta=-\frac{\pi}{r^*} \sqrt{\frac{v}{a}}, \operatorname{Pr}=\frac{v}{\alpha_m}, \Gamma_T=\lambda_E a, N b=\frac{\tau D_B}{v}\left(C_w-C_{\infty}\right), \\ & N t=\frac{\tau D_T}{v T_{\infty}}\left(T_w-T_{\infty}\right), Q=\frac{Q_0}{\rho c_p a}, \Gamma_C=\lambda_C a, R=\frac{4 T_{\infty}^3 \sigma^*}{k k^*}, L e=\frac{\alpha_m}{D_B}, B i=\frac{h_f}{k} \sqrt{\frac{v}{a}}, \varepsilon=\frac{b}{a}.\end{aligned}$

The physical quantities “skin friction”, “local Nusselt number”, and “local Sherwood number” are defined as follows:

$\begin{aligned} C_{f x} & =\frac{\tau_{w x}}{\rho u_w^2}, \quad C_{g y}=\frac{\tau_{w y}}{\rho v_w^2}, \quad N u_x=\frac{x q_w}{k\left(T_w-T_{\infty}\right)}, S h_x=\frac{x q_m}{D_B\left(C_w-C_{\infty}\right)}, \\ \tau_{w x} & =\mu\left[\frac{\partial u}{\partial z}+\frac{\Gamma}{\sqrt{2}}\left(\frac{\partial u}{\partial z}\right)^2\right]_{z=0}, \tau_{w y}=\mu\left[\frac{\partial v}{\partial z}+\frac{\Gamma}{\sqrt{2}}\left(\frac{\partial v}{\partial z}\right)^2\right]_{z=0}, q_w=\left(-\left(k+\frac{16 \sigma^* T_{\infty}^3}{3 k^*}\right)\left(\frac{\partial T}{\partial z}\right)\right)_{z=0}, \\ q_m & =-D_B\left(\frac{\partial C}{\partial z}\right)_{z=0} .\end{aligned}$

$\begin{aligned} & C_{f x} \sqrt{R e_x}=\left[f^{\prime \prime}+\frac{W}{2} f^{\prime \prime}\right]_{z=0}^2, C_{g y} \sqrt{R e_y}=\left[g^{\prime \prime}+\frac{W}{2} g^{\prime \prime}\right]_{z=0}^2, \\ & \operatorname{Re}_x^{-1 / 2} N u_x=-\left(1+\frac{4}{3} R\right) \theta^{\prime}(0), \operatorname{Re}_x^{-1 / 2} S h_x=-\phi^{\prime}(0), \\ & \operatorname{Re}_x=\frac{x u_w}{v}, \quad R e_y=\frac{y v_w}{v} .\end{aligned}$

This research examines the 3D Williamson nanofluid flow along a Riga plate in a Darcy-Forchheimer porous medium. This study utilizes the Cattaneo-Christov double-diffusion model to describe non-Fourier transport in heat and nanoparticle diffusion processes. This study also implements Buongiorno's model, which considers the impacts of Brownian motion and thermophoresis. The NDSolve function in Mathematica performs numerical evaluations of the various parameters. To confirm the precision of the adopted scheme’s, the latest outcomes of the adopted were compared with those of Loganathan et al. [47], Ragupathi et al. [49], and Umar et al. [48] in Table ~\ref{tab2}, Table ~\ref{tab3}, and Table ~\ref{tab4}. The flowchart of this approach is shown in Figure ~\ref{fig3}.

In this study, the values of factors are consistently changed for nearly all figures in order to scrutinize their impacts, while holding other factors uninterrupted throughout the analysis, such as:

$\begin{aligned} & W=H a=N t=0.2, \lambda=F r=\Gamma_T=\Gamma_C=R=Q=0.1, B i=0.3, \delta=0.4, \varepsilon=N b=0.5, \\& \operatorname{Pr}=3.0, L e=1.0 .\end{aligned}$

Table ~\ref{tab5} describes the dimensionless parameter.

Figure ~\ref{fig4}, Figure ~\ref{fig5}, Figure ~\ref{fig6}, and Figure ~\ref{fig7} show the influence of $W$ and $H a$ on the various profiles. In the presence of a Williamson fluid, the effective viscosity is expressed as a nonlinear relationship with the shear rate, where the increase in $W$ reveals a weaker shear thinning effect and an increase in the apparent viscosity in the flow domain. This high viscosity improves the inner resistance of the fluid, thereby inhibiting both axial and transverse motion and resulting in a reduction in $f^{\prime}$ and $g^{\prime}$. The diminished convective transport as a result of less fluid moment will also reduce the bearing power of the flow, in terms of thermal energy and transporting nanoparticles away. Thus, the contribution of the thermal and solutal diffusive phenomena is enhanced so that the thermal and concentration boundary layers are thickened. This implies an enhancement in the distributions of $\theta$ and $\phi$ through the boundary layer zone. It is also noted that the Riga plate increases the Lorentz force, thereby accelerating the fluid in the $\mathcal{X}$-direction, implying an improvement in $f^{\prime}$ due to the increased $H a$. Nevertheless, this greater electromagnetic forcing overrides the cross-flow in the $y$-direction owing to the redistribution of the momentum. This implies a reduction in $g^{\prime}$. For 6 and $\phi$, enhanced convection thins the boundary layers, causing faster diminishing of ambient values. The velocity distribution increases by approximately 20–30\% with a rise in the Hartmann number owing to the Lorentz force, while the corresponding temperature and concentration profiles decrease by about 10–12\% an 8–10\%, respectively. Electromagnetic forcing alters primary and secondary velocities differently because the Riga plate generates a Lorentz force primarily in the x-direction, accelerating axial flow while redistributing momentum to suppress cross-flow in the y-direction due to induced drag.

Figure ~\ref{fig8} and Figure ~\ref{fig9} show the impact of $\lambda$ and Fr on $f^{\prime}$ and $g^{\prime}$. The figures show that when $\lambda$ increases, it creates more resistance inside the porous matrix, which lowers $f^{\prime}$ and $g^{\prime}$, while making the momentum boundary layer thinner. The acceleration in Fr represents the inertial drag in the porous medium, which restricts fluid movement by adding extra resistance to the flow. The combination of porous drag and inertial resistance creates a significant velocity reduction close to the Riga plate, which shows that permeability control and inertial losses play essential roles in controlling the transport processes in non-Newtonian nanofluid systems. Physically, Darcy resistance dissipates momentum through viscous drag, while Forchheimer resistance introduces inertial losses; together, they weaken fluid motion and redistribute momentum within the boundary layer.

Figure ~\ref{fig10} and Figure ~\ref{fig11} show the repercussions of the dimensionless properties of $\mathcal{E}$ and $\delta$ on $f^{\prime}$ and $g^{\prime}$. As $\delta$ increases, the force due to electromagnetism declines very quickly away from the plate, which decreases the acceleration of the fluid along the stretching axis; as a result, $f^{\prime}$ diminishes. However, $g^{\prime}$ is enhanced owing to the relative rearrangement of the momentum in that orientation. It is also noted that the stretching ratio parameter decreases $f^{\prime}$ owing to the rise in the value bor decay of $G$. This leads to an enhancement in $g^{\prime}$.

The effects of $\Gamma_T$ and $\operatorname{Pr}$ on $\theta$ are shown in Figure ~\ref{fig12}. Enhancing parameter $\Gamma_T$ decreases the adequacy of thermal energy movement within the liquid, which implies a reduction in $\theta(\zeta)$. On the same note, an increased $\operatorname{Pr}$ signifies reduced thermal diffusivity compared with momentum diffusivity, which inhibits the transfer of heat into the fluid and the resulting decreased $\theta$.

The effects of $\varepsilon$ and $Q$ on $\theta$ are illustrated in Figure ~\ref{fig13}. A positive change in $\varepsilon$ increases the acceleration of the fluid along the surface, which reinforces the convective heat transport away from the wall, thereby decreasing $\theta$. Conversely, an elevated value of $Q$ provides an internal source of heat to the fluid, and consequently enhances the level of energy concentration in the boundary layer, resulting in an increase in $\theta$.

Figure ~\ref{fig14} show the impact of $R$ and $B i$ on $\theta$. A positive indicates $B i$ that the surface is convectively heated. In this case, a larger value of $R$ boosts the radiative heat transmission into the fluid. This extra energy amplifies the thermal boundary layer temperature, causing the temperature profile to grow. On the other hand, a negative $B i$ corresponds to convective cooling of the surface, and intensive radiation contributes to more heat transfer by the fluid. This increases the thermal energy extraction, leading to a reduced temperature profile.

Figure ~\ref{fig15} and Figure ~\ref{fig16} show the impacts of $N b$ and $N t$ on 6 and $\phi$. Th parameter $N t$ increases the migration process of nanoparticles in warmer to cooler regions thus limiting the outward heat conduction to the wall and thickens both the thermal and concentration boundary layers thus rajising both 6 and $\phi$. Contrastingly, an increased value of $N b$ leads to better diffusion of the nanoparticle in random fashion that facilitates the transfer of the thermal energy and increases 6 but at the same moment disperses more efficiently into the bulk flow, declining their near-wall concentration and thus diminishing $\phi$. The Brownian motion parameter decreases the temperature field by 8-12\% but raises the concentration field by 6-9\%, while the thermophoresis parameter enhances both temperature and concentration distributions by 15-20\% and $10-14 \%$, respectively.

Figure ~\ref{fig17} shows changes $\phi$ caused by $\Gamma_C$ and $L e$. A large $L e$ implies that the diffusivity of mass to that of temperature is very small, therefore, mass movement of the boundary layer is not large and concentrations are shallow. Similarly, the response of the nanoparticle mass flux to concentration gradients is slowed as $\Gamma_C$ is increased which dilutes the strength of diffusive transport and further causes a lowering of $\phi$.

A 2D plot can only show how one parameter changes, while keeping the other parameters same. To analyze the combined effects of two or more parameters, multiple curves or graphs are required appropriate visualization. The parameter space appears in full through the 3D surface maps that reveal essential performance regions, including peaks, valleys, ridges, and saddle points. The 3D visualization presents parameter relationships in a clearer manner, which helps users better understand the findings and conduct comparisons between results. From Figure ~\ref{fig18}, it is clear that $C_{f x}$ increases with an increase in $W$ and $H a$. Figure ~\ref{fig19} shows that $C_{g y}$ declines with both $\lambda$ and Fr. Figure ~\ref{fig20} shows that the local Nusselt number was enhanced by $H a$ and $\Gamma_T$. Figure ~\ref{fig21} illustrates that the local Sherwood number accelerates with $\Gamma_C$ and $N t$. The Brownian motion parameter slightly decreases local Nusselt number $(\approx 10–12 \%)$ but increases local Sherwood number by $7–9 \%$. Thermophoresis and Brownian diffusion play a crucial role in governing nanoparticle transport in the vicinity of the wall. Brownian diffusion promotes random motion of nanoparticles, enhancing particle mixing and increasing thermal energy transport near the surface. As a result, higher Brownian motion intensifies the heat transfer rate and elevates the local Nusselt number. Conversely, thermophoresis induces a migration of nanoparticles from the hot wall toward cooler regions, which reduces the near-wall nanoparticle concentration and weakens the mass transfer rate, leading to a decrease in the Sherwood number. The interplay between these two mechanisms controls the thickness of the thermal and concentration boundary layers and significantly influences near-wall heat and mass transfer characteristics.

Coefficients of skin friction, heat and mass transfers are highlighted in Table ~\ref{tab6}.

Since the physical boundary conditions are prescribed at infinity, the semi-infinite domain $0<\zeta<\infty$ is truncated to a finite computational domain $0<\zeta<\zeta_{\infty}$. To determine a suitable value of $\zeta_{\infty}$, numerical solutions were obtained for increasing values of $\zeta_{\infty}$ with the base set values of the parameter which were given above.

The results listed in Table ~\ref{tab7} reveal that further refinement of the computational mesh leads to insignificant changes in the values of the skin friction coefficient, Nusselt number, and Sherwood number. This confirms that the numerical solution is independent of the mesh size, and hence the adopted grid provides a stable and convergent solution.

In order to examine the effect of mesh density on the numerical solution, a mesh-independence study was conducted by varying the number of grid points (or maximum steps) used in the numerical integration.

The results listed in Table ~\ref{tab8} reveal that further refinement of the computational mesh leads to insignificant changes in the values of the skin friction coefficient, Nusselt number, and Sherwood number. This confirms that the numerical solution is independent of the mesh size, and hence the adopted grid provides a stable and convergent solution.

The response surface regression analysis was conducted to evaluate the impact of the thermal radiation parameter, Eckert number and thermal slip factor on the response variable. The model's statistical significance was validated with an $\mathrm{R}^2$ value of $99.85 \%$, an adjusted $\mathrm{R}^2$ of $99.43 \%$ and a predicted $\mathrm{R}^2$ of $99.17 \%$, demonstrating an exceptional model fit. ANOVA results indicated that all linear, quadratic and interaction terms were statistically significant, with $\mathrm{P}<0.05$. Among the linear terms, had the most significant effect $(\mathrm{F}=427.70)$, followed by $\mathrm{Ha}(\mathrm{F}=65.16)$ and $\mathrm{Q}(\mathrm{F}=2857.88)$. The quadratic term Q exhibited the highest influence $(F=294.52)$, indicating a strong nonlinear effect of. Interaction effects were substantial, with $\mathrm{R}^* \mathrm{Q}(\mathrm{F}=117.87)$ being the most dominant, followed by $\mathrm{Ha}^* \mathrm{Q}(\mathrm{F}=21.76)$ and $\mathrm{R}^* \mathrm{Ha}(\mathrm{F}=3.46)$. The final regression equation effectively predicts system behavior, with a nonsignificant lack-of fit confirming model adequacy (see Table ~\ref{tab9} and Table ~\ref{tab10}).

$N u_x=\alpha_0+\alpha_1 x+\alpha_2 y+\alpha_3 z+\alpha_{11} x^2+\alpha_{22} y^2+\alpha_{33} z^2+\alpha_{12} x y+\alpha_{13} x z+\alpha_{23} y z$;

$\begin{aligned} & N u_x=18.54+11.13 R-4.82 H a-34.86 Q-0.27 R^* R-0.43 H a^* H a+18.73 Q * Q-2.14 R^* H a \\ & -7.064 R^* Q+5.72 H a * Q;\end{aligned}$

$\begin{aligned} & N u_x=18.54+11.13 \mathrm{R}-4.82 \mathrm{Ha}-34.86 \mathrm{Q}+18.73 \mathrm{Q}^* \mathrm{Q}-2.14 \mathrm{R}^* \mathrm{Ha}-7.064 \mathrm{R}^* \mathrm{Q}+5.72 \mathrm{Ha}^* \mathrm{Q}.\end{aligned}$

The contour and surface plot Figure ~\ref{fig22}, Figure ~\ref{fig23}, and Figure ~\ref{fig24} illustrates the variation of the Nusselt number $N u$ with $H a$ and $R$. The Nusselt number ($N u$) improves with increases with respect to the Hartmann number ($H a$) and the radiation parameter ($R$), showing that greater viscous dissipation and thermal radiation improve the convective heat transfer rate in the flow.

In addition to the truncation and mesh-independence tests, residual decay curves were examined to further validate the convergence of the numerical scheme.

Figure ~\ref{fig25} illustrates the residual decay of the momentum equation obtained by substituting the numerical solution into the transformed governing equation. It is evident that the residual decreases monotonically and reaches the order of $10^{-8}$, confirming the accuracy, stability, and convergence of the numerical method employed.

The combined effects of Williamson rheology, electromagnetic forcing, porous-medium inertia, and non-Fourier fluxes are expected in engineering contexts such as MHD-assisted cooling in electronics (where $H a$ = 0.5–2.0 enhances heat removal by 8–12\%), polymer extrusion processes involving shear-thinning fluids, and biomedical transport like magnetic drug delivery in porous tissues. Parameters with the strongest design implications include the modified Hartmann number (Ha, boosting skin friction by 20–30\%) and thermophoresis ($N t$, altering concentration by 10–14\%). The present results inform optimization of cooling devices by using Riga plates to counteract porous drag, microfluidic platforms through porosity control for flow rate stability, and biomedical systems by predicting nanoparticle distribution for targeted hyperthermia or tissue engineering. The present study analyzed the 3D Williamson nanofluid flow over a stretching Riga plate embedded in a Darcy-Forchheimer porous medium, incorporating the Cattaneo-Christov heat and mass flux models along with the effects of Brownian motion, thermophoresis, and other relevant parameters. The key findings are:

• The Williamson parameter increases the effective viscosity, leading to a reduction in velocities in both orientations, while the temperature and concentration profiles rise owing to reduced convective transport and stronger diffusion.

• The modified Hartmann number enhances electromagnetic forcing, which accelerates the primary velocity but reduces the secondary velocity, temperature, and concentration distributions.

• The porous parameter and Forchheimer number act as momentum sinks, decreasing the velocities in both directions.

• The radiation parameter increases the temperature for a positive Biot number (convective heating) but decreases the temperature for a negative Biot number (convective cooling).

• The Lewis number and mass thermal relaxation time parameter reduce the concentration profiles by weakening mass diffusion within the boundary layer.

• The thermophoresis factor amplifies both $\theta$ and $\phi$, whereas the rownian motion factor amplifies $\theta$ but reduces $\phi$.

• The skin friction coefficient in x direction typically rises with stronger electromagnetic effects.

• Applying models such as Cattaneo-Christov for heat and mass flux in Williamson nanofluid flow leads to reduced thermal and concentration layer thicknesses with higher relaxation parameters, thereby improving control over heat and mass transfer rates.

• Future Directions: Further studies are recommended to explore variable viscosity, complex geometries, or hybrid nanofluid formulations to optimize the Williamson nanofluid performance in advanced engineering applications.