This article deals with developing a solution approach, called the non-isothermal negative saturation (NegSat) solution approach. The NegSat solution approach solves efficiently any non-isothermal compositional flow problem that involves phase disappearance, phase appearance, and phase transition. The advantage of the solution approach is that it circumvents using different equations for single-phase and two-phase regions and the ensuing unstable procedure. This paper shows that the NegSat solution approach can also be used for non-isothermal systems. The NegSat solution approach can be implemented efficiently in numerical simulators to tackle modeling issues for mixed CO2–water injection in geothermal reservoirs, thermal recovery processes, and for multicontact miscible and immiscible gas injection in oil reservoirs. We illustrate the approach by way of example to cold mixed CO2–water injection in a 1D geothermal reservoir. The solution is compared with an analytical solution obtained with the wave-curve method (method of characteristics) and shows excellent agreement. A complete set of simulations is carried out, which identifies six bifurcations. The two main bifurcations are (1) when the most downstream compositional wave is replaced by a compositional shock and (2) when an extra Buckley–Leverett rarefaction appears. The plot of the useful energy (exergy) versus the CO2 storage capacity shows a Z-shape. The top horizontal part represents a branch of high exergy recovery/relatively lower storage capacity, whereas the bottom part represents a branch of lower exergy recovery/higher storage capacity.