**3. Computational Methodology**

Following the HAM, choosing the initial guesses and linear operators for the velocities, temperature and homogeneous-heterogeneous chemical concentration profiles as

$$f\_0(\zeta) = \zeta^3(k\_6 + k\_8) - \zeta^2(2k\_6 + k\_8) + \zeta k\_6, \quad g\_0(\zeta) = \zeta \Omega + 1 - \zeta, \quad \theta\_0(\zeta) = -\zeta + 1,$$

$$\rho\_0(\zeta) = \frac{\zeta k\_7 + 1}{k\_7 + 1}, \quad \text{(41)}$$

$$\mathbf{g}^{\prime\prime} = \mathbf{L}\_{\mathbf{\varphi}\prime} \qquad \qquad \mathbf{f}^{\prime\prime\prime} = \mathbf{L}\_{\mathbf{f}\prime} \qquad \qquad \mathbf{g}^{\prime\prime} = \mathbf{L}\_{\mathbf{\mathcal{G}}\prime} \qquad \qquad \boldsymbol{\theta}^{\prime\prime} = \mathbf{L}\_{\boldsymbol{\theta}\prime} \tag{42}$$

characterizing

$$\begin{array}{rclcrcl} \mathbf{L}\_f \begin{bmatrix} \mathbf{E}\_1 + \mathbf{E}\_2 \boldsymbol{\zeta} + \mathbf{E}\_3 \boldsymbol{\zeta}^2 + \mathbf{E}\_4 \boldsymbol{\zeta}^3 \end{bmatrix} &=& \mathbf{0}, & \mathbf{L}\_\mathbf{\boldsymbol{\zeta}} \begin{bmatrix} \mathbf{E}\_5 + \mathbf{E}\_6 \boldsymbol{\zeta} \end{bmatrix} &=& \mathbf{0}, & \mathbf{L}\_\mathbf{\boldsymbol{\theta}} \begin{bmatrix} \mathbf{E}\_7 + \mathbf{E}\_8 \boldsymbol{\zeta} \end{bmatrix} &=& \mathbf{0}, & \mathbf{L}\_\mathbf{\boldsymbol{\theta}} \begin{bmatrix} \mathbf{E}\_9 + \mathbf{E}\_{10} \boldsymbol{\zeta} \end{bmatrix} &=& \mathbf{0}, & \mathbf{\{43\}} \end{array}$$

where *E<sup>i</sup>* (*i* = 1–10) are the arbitrary constants.
