A boundary element method has been developed for analysing heat transport phenomena in solitarywave on falling thin liquid films at high Reynolds numbers. The divergence theorem is applied to thenon-linear convective volume integral of the boundary element formulation with the pressure penaltyfunction. Consequently, velocity and temperature gradients are eliminated, and the complete formulationis written in terms of velocity and temperature. This provides considerable reduction in storage andcomputational requirements while improving accuracy. The non-linear equation systems of boundaryelement discretization are solved by the quasi-Newton iterative scheme with Broyden’s update. Thestreamline maps and the temperature distributions in solitary wave and wavy film flow have been obtained,and the variations of Nusselt numbers along the wall-liquid interface are also given. There arelarge cross-flow velocities and S-shape temperature distributions in the recirculating region of solitarywave. This special flow and thermal process can be a mechanism to enhance heat transport. A boundary element method has been developed for analysing heat transport phenomena in solitarywave on falling thin liquid films at high Reynolds numbers. The divergence theorem is applied to then-linear convective volume integral of the boundary element formulation with the pressure penalty function. This non-linear equation systems of boundaryelement discretization are solved by the quasi-Newton iterative scheme with Broyden’s update There are large cross-flow velocities and S-shape temperature distributions in the recirculating region. of solitarywave. This special flow and thermal process can be a mechanism to enhance heat transport.