We report results of long timescale adaptive kinetic Monte Carlo simulations aimed at identifying possible molecular reordering processes on both proton-disordered and ordered (Fletcher) basal plane (0001) surfaces of hexagonal ice. The simulations are based on a force field for flexible molecules and span a time interval of up to 50 μs at a temperature of 100 K, which represents a lower bound to the temperature range of earth's atmosphere. Additional calculations using both density functional theory and an ab initio based polarizable potential function are performed to test and refine the force field predictions. Several distinct processes are found to occur readily even at this low temperature, including concerted reorientation (flipping) of neighboring surface molecules, which changes the pattern of dangling H-atoms, and the formation of interstitial defects by the downwards motion of upper-bilayer molecules. On the proton-disordered surface, one major surface roughening process is observed that significantly disrupts the crystalline structure. Despite much longer simulation time, such roughening processes are not observed on the highly ordered Fletcher surface which is energetically more stable because of smaller repulsive interaction between neighboring dangling H-atoms. However, a more localized process takes place on the Fletcher surface involving a surface molecule transiently leaving its lattice site. The flipping process provides a facile pathway of increasing proton-order and stabilizing the surface, supporting a predominantly Fletcher-like ordering of low-temperature ice surfaces. Our simulations also show that eventual proton-disordered patches on the surface may induce significant local reconstructions. Further, a subset of the molecules on the Fletcher surface are susceptible to forming interstitial defects which might provide active sites for various chemical reactions in the atmosphere.