Nanoscale multilayers of binary metallic systems, such as nickel/aluminum, exhibit self-propagating exothermic reactions due to the high formation enthalpy of the intermetallic compounds. Most of the previous modeling approaches on the reactions of this system rely on the use of mass diffusion with a phenomenological derived diffusion coefficient representing single-phase (NiAl) growth, coupled with heat transport. We show that the reaction kinetics, temperatures, and thermal front width can be reproduced more satisfactorily with the sequential growth of Ni2Al3 followed by NiAl, utilizing independently obtained interdiffusivities. The computational domain was meshed with a dynamically generated bi-modal grid consisting of fine and coarse zones corresponding to rapid and slower reacting regions to improve computational efficiency. The PDEPE function in MATLAB was used as a basis for an alternating direction scheme. A modified parabolic growth law was employed to model intermetallic growth in the thickness direction. A multiphase enthalpy function was formulated to solve for temperatures after discrete phase growth and transformations at each time step. The results show that the Ni2Al3 formation yields a preheating zone to facilitate the slower growth of NiAl. At bilayer thicknesses lower than 12 nm, the intermixing layer induces oscillating thermal fronts, sharply reducing the average velocities. (C) 2015 AIP Publishing LLC.