The shallow-ice approximation is used to analyse the normal-mode stability of long-wavelength, thermally coupled, free-surface flows of glaciers and ice sheets along infinite flow sections. The viscosity of the ice changes with temperature, and the ice is also thermally coupled with the bedrock. An assumption of quasi-uniform flow, where the downstream flux of mass and energy is assumed to be constant over relatively short wavelengths, permits vertical advection due to accumulation to be considered. The assumption allows realistic representation of the vertical velocity distribution and temperature distributions within the modelled flow. Nevertheless, the linearized equations are formally independent of the assumption of quasi-uniform flow. The linearized equations are Fourier-transformed over the horizontal plane and an algebraic eigenvalue problem constructed which considers the Fourier-transforms of the vertically varying temperature and thickness. The dependence of the eigenvalues on wavenumber, both parallel and orthogonal to flow, accumulation rate, geothermal heat flux, surface temperature and slope is computed. Extensive linearly unstable regimes where the base is well below melting point are found, which correspond only approximately to negative slopes in the flux/thickness relationship (supercritical zones). The most unstable parts of these regimes have finite along-flow wavelength, and have very long transverse wavelength. These modes are found to be more unstable in two senses: (i) they bifurcate for lesser thicknesses; (ii) they have greater growth rates. All computed wave velocities of unstable modes were downslope. Less unstable modes with short and medium transverse wavelengths are found which excite the temperature field only. Time constants range between ten and one thousand years. There is no strong evidence for stream-pattern formation through a thermoviscous instability mechanism.