Numerical simulation of Taylor flows presents several challenges. At the dynamic interface physical properties are discontinuous, which is especially challenging for the thin film between the droplet and the wall. Phase-field methods, which are derived from thermodynamic principles, define the interface as a smooth transition between phases. By coupling the Cahn-Hilliard equation with the Navier-Stokes and energy equation, both interface dynamics and heat transfer can be captured. In the work presented, the resulting system of equations are solved by a parallel h-adaptive least-squares spectral element method. To approximate the solution with sufficient numerical accuracy, C1 Hermite basis functions and a space-time formulation have been applied. It is widely accepted in the literature that the droplet characteristics such as length, velocity and dynamic interaction among them affect the heat transfer properties of Taylor flow. To gain understanding, their effect on heat transfer and pressure drop for liquid-liquid Taylor flow in microchannels must be studied in more detail.