Charge trapping and formation of polarons is a pervasive phenomenon in transition metal oxide compounds, in particular at the surface, affecting fundamental physical properties and functionalities of the hosting materials. Here we investigate via first-principle techniques the formation and dynamics of small polarons on the reduced surface of titanium dioxide, an archetypal system for polarons, for a wide range of oxygen vacancy concentrations. We report how the essential features of polarons can be adequately accounted in terms of few quantities: the local structural and chemical environment, the attractive interaction between negatively charged polarons and positively charged oxygen vacancies, and the spin-dependent polaron-polaron Coulomb repulsion. We combined molecular dynamics simulations on realistic samples derived from experimental observations with simplified static models, aiming to disentangle the various variables at play. We find that depending on the specific trapping site, different types of polarons can be formed, with distinct orbital symmetries and different degree of localization and structural distortion. The energetically most stable polaron site is the subsurface Ti atom below the undercoordinated surface Ti atom, owing to a small energy cost to distort the lattice and a suitable electrostatic potential. Polaron-polaron repulsion and polaron-vacancy attraction determine the spatial distribution of polarons as well as the energy of the polaronic in-gap state. In the range of experimentally reachable oxygen vacancy concentrations the calculated data are in excellent agreement with observations, thus validating the overall interpretation.arXiv:1805.01849v1 [cond-mat.mtrl-sci]