By solving a 1D heat conservation equation for single phase flow, Butler (1981) derived his classical SAGD equation which has excellent predictive capability at experimental scales but performs poorly at field scales. Several authors have tried to remedy this by acknowledging the existence of multiphase flow ahead of the steam chamber during SAGD, but have only empirically accounted for it using multipliers that seem to vary for each reservoir. Recently, by coupling the mass and energy conservation equations, Azom et al. (2013) solved this problem and showed that the multi-scale phenomenon is controlled by the classical Marangoni or thermo-capillary number. In this paper, we extend the concept of thermo-capillary imbibition (Azom et al., 2013) to include the additional coupling of a solvent component mass balance to the mass and energy conservation equations in order to derive an expression for oil recovery during ES-SAGD. Our results show that in addition to the thermo-capillary number, additional dimensionless groups like the condensate Lewis number, the ratios of the solvent's condensate to bitumen dispersion coefficient and solvent viscosity to bitumen viscosity at steam temperature influence multiphase flow and hence ES-SAGD flow performance. An important conclusion is that heterogeneity will have a complex effect on the ES-SAGD process as it causes more thermo-capillary imbibition (thereby decreasing bitumen rates) while at the same time causes the further transport of solvent into the bitumen phase (thereby increasing bitumen rates). Our method involved accounting for capillary pressure in the conservation equations and decoupling it into its temperature and saturation components. This then forms the basis for coupling the mass and energy conservation equations. The solvent is assumed not to affect capillary behavior for the ES-SAGD process. The entire modeling was done in dimensionless space making the results applicable across varying scales and for a range of reservoir parameters. This work presents a viable semi-analytical model for the ES-SAGD process while at the same time accounting for the effect of multiphase flow during the process. It is also important to note that currently no commercial simulator can model thermocapillary behavior and hence cannot account for multiphase flow effects during non-isothermal recovery processes. Our model solves this problem for the ES-SAGD process and can be used together with the original Butler model as a fast ES-SAGD predictive model that also accounts for multiphase flow in any proxy-based history matching process.