<p>A wide spectrum of natural phenomena and industrial processes involve the propagation of reaction fronts in a viscous fluid medium. Consistent models of this class of problems should account for the competition between the flow dynamics and the reaction properties related to temperature and degree of conversion. In the present study, to resolve these features in reactive flows, we consider the incompressible Navier–Stokes equations with temperature-dependent viscosity for the velocity and pressure coupled with a double system of reaction-advection-diffusion equations for the temperature and the degree of conversion. The reactions in this coupled system are assumed to be generated using the Arrhenius activation energy and cross diffusion coefficients are influenced by both the temperature and degree of conversion. We analyze the solvability of the proposed model and prove existence of the weak solution using a Galerkin approach combined with appropriate compactness arguments. The uniqueness of the solution is also achieved under additional assumptions on the problem data. In the current work, we also introduce a time-marching scheme for the numerical approximation of the continuous problem and prove the unconditional stability with respect to the time step. Numerical results are presented for a problem with exact solutions to verify the theoretical analysis and to assess the performance of the numerical method. The coupled system is also employed to solve a two-dimensional flame-like propagation problem in viscous fluids. The computational results obtained for both examples support the theoretical expectations for a stable and accurate numerical solver for reactive fluids with cross diffusion and Arrhenius activation energy.</p>