This work introduces a highly accurate and efficient splitting method for solving thermally coupled incompressible flow problems involving nonlinear viscoelastic fluids. The proposed technique employs a semi-implicit reformulation of the governing system, enabling the efficient decoupling of individual sub-problems (velocity, pressure, stresses, and temperature) at each time step. A notable strength of the framework is its high-order accuracy and its flexibility to accommodate various viscoelastic constitutive models [1]. The method is evaluated using the benchmark problem of a differentially heated square cavity, which incorporates viscous dissipation and an enthalpy-based phase-change model [2]. Numerical experiments are conducted for Phan-Thien-Tanner (PTT) and FENE-P fluids, combined for the first time with a phase-change model. The study considers Rayleigh numbers in the range of 10^3 to 10^5 and Weissenberg numbers up to 10, exploring the effects of flow parameters on the hydrodynamic and thermal boundary layer patterns. Average Nusselt number values are validated against reference results, while additional findings are presented to guide future research. Qualitative results for velocity, temperature, and stress fields are also included, demonstrating the method’s capability to capture intricate details of natural convection and phase-change cycles under both steady and transient conditions. The results provide a comprehensive description of fluid flow and the solid-liquid-solid phase-change process within the cavity, offering valuable insights into these complex thermal and mechanical phenomena.