## 3 Solution of the Problem  
The boundary value problem given in Equations (19)-(25) admits the closed form analytical solutions for the temperature, concentration, and velocity profiles in the following way  
Using Equations (22) and the boundary conditions for \(\theta\) defined in Equations (24) and (25), we obtain the following expression for the temperature distribution:  
\[\theta=1-\frac{S_{t}}{2}\Big{(}y^{2}-h^{2}\Big{)}. \tag{26}\]  
Differentiating Equation (26) twice with respect to \(y\) and using the result in Equation (23), we obtain:  
\[\frac{1}{S_{c}}\frac{\partial^{2}\varphi}{\partial y^{2}}-M\varphi-S_{r}S_{t }=0, \tag{27}\]  
which leads to the following solution satisfying the relevant boundary for \(\varphi\):  
\[\varphi=\frac{\cosh\Big{(}\sqrt{MS_{c}}y\Big{)}}{\cosh\Big{(}\sqrt{MS_{c}}h \Big{)}}\Big{(}1+\frac{S_{r}S_{t}}{M}\Big{)}-\frac{S_{r}S_{t}}{M}. \tag{28}\]  
Upon using Equations (26) and (28) and boundary conditions, the axial velocity distribution \(u\) from Equation (20) is formulated as:  
\[u=\Big{(}-1-\frac{2\pi\varepsilon\delta\cos(2\pi x)}{1-2\pi\varepsilon\delta \cos(2\pi x)}\Big{)}\frac{\cosh\Big{(}\sqrt{MS_{1}}y\Big{)}}{\cosh\Big{(} \sqrt{MS_{1}}x\Big{)}}+X_{3}\Big{(}\frac{\cosh\Big{(}\sqrt{MS_{1}}y\Big{)}}{ \cosh\Big{(}\sqrt{MS_{1}}x\Big{)}}-\frac{\cosh\big{(}\sqrt{MS_{c}}y\big{)}}{ \cosh\Big{(}\sqrt{MS_{1}}x\Big{)}}\Big{)}+\Big{[}X_{4}+X_{5}\frac{d\varphi}{ \varepsilon\varepsilon\varepsilon}\Big{(}\frac{\cosh\big{(}\sqrt{MS_{1}}y \big{)}}{\cosh\Big{(}\sqrt{MS_{1}}x\big{)}}-1\Big{)}+X_{6}\Big{(}y^{2}-h^{2} \Big{)} \tag{29}\]  
where \(X_{1}=(1+\lambda_{1})\Big{(}H_{a}^{2}+\frac{1}{D_{a}}\Big{)}\), \(X_{2}=X_{1}+(1+\lambda_{1})\frac{\delta\varepsilon S_{t}S_{t}}{M}\), \(X_{3}=\frac{(1+\lambda_{1})S_{c}}{(MS_{c}-X_{1})}\Big{(}1+\frac{S_{r}S_{t}}{M} \Big{)}\), \(X_{4}=\frac{X_{2}+\delta S_{t}S_{t}(1+\lambda_{1})}{X_{1}}\), \(X_{5}=\frac{(1+\lambda_{1})}{(MS_{c}-X_{1})}\), and \(X_{6}=\frac{-(1+\lambda_{1})S_{t}S_{t}}{2X_{1}}\).  
In order to determine the unknown pressure gradient in Equation (28), we used the following relation:  
\[q=\int\limits_{0}^{h}udy, \tag{30}\]which upon utilizing Equation (29) and then solving for \(\frac{dp}{dx}\), takes the following form:  
\[\frac{dp}{dx}=\frac{1}{\tanh\sqrt{X_{1}}h-\sqrt{X_{1}}h}\Bigg{[}\frac{\sqrt{X_{1}} q}{X_{5}}-\frac{X_{6}{\tanh\sqrt{X_{1}}}h}{X_{5}}-\frac{X_{3}{\tanh\sqrt{X_{1}}}h}{X_{5}}+ \frac{X_{3}\sqrt{X_{1}}{\tanh\sqrt{MS_{c}}}h}{X_{5}\sqrt{MS_{c}}}\Bigg{]}, \tag{31}\]  
where \(q\) is the volume flow rate, which is interconnected with \(Q\) (the nondimensional time mean volume flow rate) as follows:  
\[q=Q-1. \tag{32}\]  
Upon integrating Equation (31), we obtain \(\Delta p\) (the pressure rise or drop) over one wavelength as:  
\[\Delta p=\int\limits_{0}^{1}\frac{dp}{dx}dx=a_{11}(Q-1)I_{1}-a_{12}I_{2}-a_{13  
parameter). Figure 3c explains the velocity rises when we increased the value of the parameter \(\varepsilon\). The overall features of the velocity profile in the present study for all cases were in good agreement with the existing studies in the literature (17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28); for instance, the flow was symmetric, the velocity profiles exhibited parabolic graphs, their maximum values were almost near the central region and decreased as we moved from center to the walls of the channel. This validates the correctness and generality of our model.