The thermal evolution of Ganymede and Callisto was calculated using a realistic water ice heat flux scaling law. This scaling was derived from numerical studies of thermal convection with a stress and temperature dependent viscosity relationship that is descriptive of dislocation creep of crystalline water ice I. This result is then compared with that determined using a non-Newtonian silicate like rheology, a strongly temperature dependent Newtonian rheology and an isoviscous scaling law. In all cases, the stagnant lid convection models resulted in warmer internal mantle temperatures and thicker conductive lids than the isoviscous scaling law, in agreement with theory. It was found that the isoviscous thermal evolution models for each icy satellite gave cool, frozen models, whilst the stagnant lid scaling laws resulted in warmer satellite interiors. High viscosity models resulted in extensive melting of the icy mantle. The presence of subsurface oceans within the icy Galilean satellites has been inferred from the Galileo data and stagnant lid thermal convection is one mechanism through which a subsurface ocean may form.