Graf incorporated heat transport within the saturated-zone flow regime into HydroGeoSphere together with temperature-dependent fluid properties, such as viscosity and density. The model’s capability was successfully demonstrated for the case of thermohaline flow and transport in porous and fractured porous media [Graf and Therrien, 2007. This work extends the model’s capability to include thermal energy transport in the unsaturated zone and in the surface water, which is considered a key step in the linkage between the atmospheric and hydrologic systems. Surface heat fluxes from atmospheric inputs are an important source/sink of thermal energy, especially to the surface water system. As such, surface heat fluxes across the land surface were also incorporated into HydroGeoSphere. A complete description of the physical processes and governing flow and solute transport equations that form the basis of HydroGeoSphere can be found in Therrien et al. and therefore will not be presented here.