The integral method of Thwaites, for computing the primary parameters of laminar boundary layers with constant fluid properties, is extended to heated boundary layers in water, taking into account variable fluid properties. Universal parameters are correlated from numerical solutions of heated water wedge flows for use with the integral method. The method shows good accuracy in a test with the Howarth retarded flow. The Lighthill high Prandtl number approximation is extended to permit computation of the Nusselt number for boundary layers with variable fluid properties. Nusselt numbers computed for the Howarth flow are close to the exact numerical solutions, except near separation.