In this work we propose a novel approach to model order reduction for incompressible fluid flows that focuses on the spatio-temporal description of the stresses on the surface of a body, i.e. of the wall shear stress and of the wall pressure. The spatial representation of these two variables is given by a compact set of "wall basis functions", i.e. elementary basis functions defined on the wall. In this paper, these are derived using the well-known Proper Orthogonal Decomposition, to represent optimally the fluctuation energy of the pressure and shear stress. On the other hand, the functional structure of the dynamic model is derived from first principles using the vorticity form of the Navier-Stokes equations, yielding a set of nonlinear ordinary differential equations for the time-varying amplitudes of the wall shear stress basis functions. Coefficients of this model are then identified from simulation data. To complete the system, we show that the surface pressure distribution, i.e. the time-varying amplitudes of the wall pressure basis functions, can be derived from a quadratic model of the wall shear stress temporal coefficients, stemming from the Poisson equation for the pressure. This further step is crucial for the correct representation of the aerodynamic forces. As a paradigmatic example, we present our approach for the modelling of the free dynamics of the separated flow around a circular cylinder in the laminar regime, at Re = 200. Further implications and potentialities of the proposed approach are discussed.