  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 {\mu}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, but our simulations
also show that 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.
