An accurate time integration method for diffusion and kinematic overland flow approximations are proposed. The discretization of first- and second-order spatial derivatives, in the basic equation is obtained by using second-order Lax-Wendroff and three-point centred finite difference schemes, respectively. For the solution in time, the system of ordinary differential equations, obtained by linearization of the celerity, of the hydraulic diffusivity and by Taylor series expansions, is integrated analytically. The stability and the numerical dissipation and dispersion are investigated by Fourier linear analysis. A proper Courant number, and the corresponding time step for the numerical simulations can be established. In addition, the proposed diffusion and kinematic wave models are straightforwardly extended to the two-dimensional flow. Test cases with both one- and two-dimensional problems, compared with analytical solutions, experimental data and with the numerical solution of full Saint-Venant equations are presented. These simulations show that the proposed numerical-analytical models accurately predict the overland flow for several situations, in particular for unsteady rainfall rate and for spatial variations in surface roughness.