In order to study surface disordering of ice at temperatures below the bulk melting point as a function of depth into the bulk, Car-Parrinello molecular-dynamics simulations of a periodic model of the hexagonal ice (0001) surface were carried out. Partial disorder in the uppermost bilayer was observed at a simulation temperature of 190 K, which is similar to 30 K below the estimated bulk melting point, qualitatively validating earlier classical molecular-dynamics studies of this phenomenon. Over 0.5 ps, the time scale of a simulation, there were three particularly useful (and complementary) measures of disorder: The pair distribution function g(r), the distance of the oxygen atoms from the bottommost bilayer, and the distribution of angles phi and theta formed by the molecular dipole vector and the Cartesian axes. Our results set the stage for future studies addressing the effect of the disordered ice surface on heterogeneous atmospheric chemistry.