We present a theoretical study of the dynamics of H atoms adsorbed on graphene bilayers with Bernal stacking. First, through extensive density functional theory calculations, including van der Waals interactions, we obtain the activation barriers involved in the desorption and migration processes of a single H atom. These barriers, along with attempt rates and the energetics of H pairs, are used as input parameters in kinetic Monte Carlo simulations to study the time evolution of an initial random distribution of adsorbed H atoms. The simulations reveal that, at room temperature, H atoms occupy only one sublattice before they completely desorb or form clusters. This sublattice selectivity in the distribution of H atoms may last for sufficiently long periods of time upon lowering the temperature down to 0 C. The final fate of the H atoms, namely, desorption or cluster formation, depends on the actual relative values of the activation barriers which can be tuned by doping. In some cases a sublattice selectivity can be obtained for periods of time experimentally relevant even at room temperature. This result shows the possibility for observation and applications of the ferromagnetic state associated with such distribution.